Final  Report  On: 


The  Geometric  Conservation  Based  Algorithms 
for  Multi-Fluid  Flow  at  All  Speeds 

(F61775-01-W-E005,  SPC  01-4005) 

Submitted  to: 

European  Office  of  Aerospace  Research  and  Development 

(EOARD) 

by 

F.  Moukalled 

American  University  of  Beirut, 

Faculty  of  Engineering  &  Architecture, 

Mechanical  Engineering  Department, 

P.O.Box  11-0236 
Beirut  -  Lebanon 

Address  all  correspondence  to  Dr.  F.  Moukalled, 
email:  memouk@aub.edu.lb 


September  10,  2002 


REPORT  DOCUMENTATION  PAGE 


Form  Approved  OMB  No.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  the  burden,  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704-0188),  1215  Jefferson 
Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply 
with  a  collection  of  information  if  it  does  not  display  a  currently  valid  OMB  control  number. 


PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 


1.  REPORT  DATE  (DD-MM-YYYY)  2.  REPORT  TYPE 

18-09-2002  Final  Report 

3.  DATES  COVERED  (From  -  To) 

14  August  2001  -  14-Aug-02 

4.  TITLE  AND  SUBTITLE 

The  Geometric  Conservation  Based  Algorithms  For  Multi-Fluid  Flow  At  All 

Speeds 

5a.  CONTRACT  NUMBER 

F61 775-01 -WE005 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

Professor  Fadl  Moukalled 

5d.  PROJECT  NUMBER 

5d.  TASK  NUMBER 

5e.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

American  University  of  Beirut 

ME  Dept.  PO  Box:  11-0236 

Beirut 

Lebanon 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

N/A 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

EOARD 

PSC  802  BOX  14 

FPO  09499-0014 

10.  SPONSOR/MONITOR’S  ACRONYM(S) 

11.  SPONSOR/MONITOR’S  REPORT  NUMBER(S) 

SPC  01-4005 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  is  unlimited. 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

This  report  results  from  a  contract  tasking  American  University  of  Beirut  as  follows:  The  proposed  work  involves  implementing  and 

(thoroughly)  testing  the  Geometric  Conservation  Based  family  of  Algorithms  within  a  structured  finite-volume  framework.  The  convection  terms 
along  the  control  volume  faces  will  be  evaluated  using  a  High  Resolution  scheme  applied  within  the  context  of  the  Normalized  Variable 
Formulation  methodology.  In  order  to  accelerate  the  convergence  rate  and  reduce  the  overall  computational  cost  of  the  algorithm,  the  outer 
iterations  will  be  accelerated  by  using  a  non-linear  full  multigrid  method.  The  discretization  scheme  will  be  second-order  accurate  in  space  and 
first  order  accurate  in  time  (even  though  extension  to  second  order  accuracy  should  be  straightforward).  The  newly  implemented  algorithms 
will  be  tested  by  solving  a  variety  of  two-dimensional  multi-fluid  flow  problems  in  the  subsonic,  transonic,  and  supersonic  regimes.  Examples 
include:  1)  Phase  separation  in  a  duct  (water-air);  2)  Turbulent  bubbly  flow  in  a  pipe;  3)  Turbulent  gas-solid  flow  in  a  curved  duct;  4)  Dusty  flow 
over  a  flat  plate  at  subsonic  flow  conditions  (-incompressible);  5)  Dusty  flow  in  a  converging-diverging  nozzle,  and  6)  Any  problem  of  interest 
to  Dr.  Sekar  (AFRL/PR).  A  detailed  report  describing  the  work  will  be  submitted  at  end  of  contract  as  specified  in  the  schedule  of  supplies  in 
the  submitted  proposal. 


15.  SUBJECT  TERMS 

EOARD,  Computational  Fluid  Dynamics  (CFD),  High  Speed  Aerodynamic,  Two-Phase  Flows 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 
ABSTRACT 

UL 

18,  NUMBER 
OF  PAGES 

19a.  NAME  OF  RESPONSIBLE  PERSON 

Wayne  Donaldson 

a.  REPORT 

UNCLAS 

b.  ABSTRACT 

UNCLAS 

c.  THIS  PAGE 

UNCLAS 

19b.  TELEPHONE  NUMBER  (Include  area  code) 

+44  (0)20  7514  4299 

Standard  Form  298  (Rev.  8/98) 


Prescribed  by  ANSI  Std.  Z39-18 


me  Lreomemc  conservation  tsasea  Aigoritnms  jor  iviuiti-t  luia  now  at  All  speeas 


z. 


Abstract 

This  work  is  concerned  with  the  implementation  and  testing,  within  a  structured  collocated 
finite-volume  framework,  of  seven  segregated  algorithms  for  the  prediction  of  multi-phase  flow 
at  all  speeds.  These  algorithms  belong  to  the  Geometric  Conservation  Based  Algorithms 
(GCBA)  group  in  which  the  pressure  correction  equation  is  derived  from  the  constraint  equation 
on  volume  fractions  (i.e.  sum  of  volume  fractions  equals  1).  The  pressure  correction  schemes  in 
these  algorithms  are  based  on  SIMPLE,  SIMPLEC,  SIMPLEX,  SIMPLEM,  SIMPLEST,  PISO, 
and  PRIME.  Solving  a  variety  of  one-  and  two-dimensional  laminar  and  turbulent  two-phase 
flow  problems  in  the  subsonic,  transonic,  and  supersonic  regimes  and  comparing  results  with 
published  numerical  and/or  experimental  data  assess  the  performance  and  accuracy  of  these 
algorithms.  The  SG  method  is  used  to  solve  for  the  one-dimensional  test  problems  and  the 
effects  of  grid  size  on  convergence  characteristics  are  analyzed.  On  the  other  hand,  solutions  for 
the  two-dimensional  problems  are  generated  for  several  grid  systems  using  the  single  grid 
method  (SG),  the  prolongation  grid  method  (PG),  and  the  full  non-linear  multi-grid  method 
(FMG)  and  their  effects  on  convergence  behavior  are  studied.  The  main  outcomes  of  this  study 
are  the  clear  demonstrations  of:  (i)  the  capability  of  all  GCBA  algorithms  to  deal  with  multi¬ 
fluid  flow  situations;  (ii)  the  ability  of  the  FMG  method  to  tackle  the  added  non-linearity  of 
multi-fluid  flows;  (iii)  and  the  capacity  of  the  GCBA  algorithms  to  predict  multi-fluid  flow  at  all 


speeds. 
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Nomenclature 

A)1,' ,..  coefficients  in  the  discretized  equation  for  (j)a). 

Bp  '  source  term  in  the  discretized  equation  for  §<k> . 

Ba  1  body  force  per  unit  volume  of  fluid/phase  k. 

Cpk>  coefficient  equals  to  1  /R{k)T(k). 

D(r,k>[<i)<k']  the  1)  operator. 

//p[<() u  '  |  the  H  operator. 

HPp[§a>]  the  HP  operator  working  on  (f)<k)  (4>fk)=ulk),vfk),  or  w(k)). 
HPP[ua  ]  the  vector  form  of  the  HP  operator. 

I(A  >  inter-phase  momentum  transfer. 

,]{k)D  diffusion  flux  of  <|)a)  across  cell  face  ‘f . 

j‘AK  convection  flux  of  across  cell  face  T. 

Ma<  <  mass  source  per  unit  volume. 

P  pressure. 

Pr<k^Prt(k,  laminar  and  turbulent  Prandtl  number  for  fluid/phase  k. 
q(k)  heat  generated  per  unit  volume  of  fluid/phase  k. 

Q[k)  general  source  term  of  fluid/phase  k. 

ra)  volume  fraction  of  fluid/phase  k. 

Ra '  gas  constant  for  fluid/phase  k. 

Sf  surface  vector, 

t  time. 

T{k  ’  temperature  of  fluid/phase  k. 

U(k>  interface  flux  velocity  (v^.S^)  of  lluid/phase  k. 

uU)  velocity  vector  of  fluid/phase  k. 

u(k),v(k),..  velocity  components  of  fluid/phase  k. 
x,  y  Cartesian  coordinates. 

||a,b||  the  maximum  of  a  and  b. 


me  Lreomemc  conservation  tsasea  Aigoritnms  jor  iviuiti-t  luia  now  at  All  speeas 


Greek  Symbols 

pa  ’  density  of  fluid/phase  k. 

ra  ’  diffusion  coefficient  of  fluid/phase  k. 

<f> <k  1  dissipation  term  in  energy  equation  of  fluid/phase  k. 

<f)a  i  general  scalar  quantity  associated  with  fluid/phase  k. 

K  f  space  vector  equal  to  (nf  -  ydf  )sf 

A  p  [(f) ,k  ’  J  the  A  operator. 

|l {k  ’ ,  pfa '  laminar  and  turbulent  viscosity  of  fluid/phase  k. 

£1  cell  volume. 

(3  a 0  thermal  expansion  coefficient  for  phase/fluid  k. 

8t  time  step. 

Subscripts 

e,  w, .  refers  to  the  east,  west,  . . .  face  of  a  control  volume. 

E,W,..  refers  to  the  East,  West,  ...  neighbors  of  the  main  grid  point, 

f  refers  to  control  volume  face  f. 

P  refers  to  the  P  grid  point. 

Superscripts 

C  refers  to  convection  contribution. 

D  refers  to  diffusion  contribution. 

(k)  refers  to  fluid/phase  k. 

( k )  *  refers  to  updated  value  at  the  current  iteration. 

(k)  o  refers  to  values  of  fluid/phase  k  from  the  previous  iteration. 

ay  refers  to  correction  field  of  phase/fluid  k. 

old  refers  to  values  from  the  previous  time  step. 


Introduction 


The  last  two  decades  have  witnessed  a  substantial  transformation  in  the  CFD  industry;  from  a 
research  means  confined  to  research  laboratories,  CFD  has  emerged  as  an  every  day 
engineering  tool  for  a  wide  range  of  industries  (Aeronautics,  Automobile,  HVAC,  etc...). 
This  increasing  dependence  on  CFD  is  due  to  a  multitude  of  factors  that  have  rendered 
practical  the  simulation  of  complex  problems.  Some  of  these  factors  are  directly  related  to 
the  maturity  of  several  numerical  aspects  at  the  core  of  CFD.  These  include:  multi-grid 
acceleration  techniques  [1-4]  with  enhanced  equation  solvers  [5,6]  that  have  decreased  the 
computational  cost  of  tackling  large  problems,  better  discretization  techniques,  unstructured 
grids  [7-12],  bounded  high  resolution  schemes  [13-18],  as  well  as  improved  pressure- velocity 
(and  density)  coupling  algorithms  for  fluid  flow  at  all  speeds  [19-27].  Other  factors, 
independent  of  the  CFD  industry,  have  to  do  with  the  exponential  increase  in  processor 
power  and  decrease  in  microprocessor  cost,  whereby  multiprocessors  systems  with  large 
memory  can  now  be  set  up  at  a  fraction  the  cost  of  the  super-computers  of  a  decade  ago. 
Challenges  still  abound  in  relation  to  increasing  the  robustness  of  numerical  techniques, 
improving  the  models  used  (e.g.  turbulence),  and  extending  the  currently  used  algorithms 
[28-34]  for  the  simulation  of  multi-phase  flows  at  all  speeds  [35].  In  this  last  area  a  number 
of  algorithms  have  been  recently  reviewed  and  new  ones  proposed  [36] .  The  basic  difficulty 
in  the  simulation  of  multi-phase  flows  [36]  stems  from  the  increased  algorithmic  complexity 
that  need  to  be  addressed  when  dealing  with  multiple  sets  of  continuity  and  momentum 
equations  that  are  inter-coupled  (interchange  momentum  by  inter-phase  mass  and  momentum 
transfer,  etc.)  both  spatially  and  across  fluids.  Despite  these  complexities,  successful 
segregated  incompressible  pressure-based  solution  algorithms  have  been  devised.  The  IPSA 
variants  devised  by  the  Spalding  Group  at  Imperial  College  [37-39]  and  the  set  of  algorithms 
devised  by  the  Los  Alamos  Scientific  Laboratory  (LASL)  group  [40-42]  are  examples  of 
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incompressible  multiphase  algorithms.  When  dealing  with  all-speed  flows,  pressure- 
velocity-density  coupling  has  to  be  accounted  for.  Pressure -based  algorithms  have  been 
extended  successfully  [19-27]  to  account  for  this  additional  coupling. 

Recently,  Darwish  et  al.  [36]  extended  the  applicability  of  the  available  segregated  single¬ 
fluid  flow  algorithms  [35]  to  predict  multi-fluid  flow  at  all  speeds.  In  their  work,  it  was 
shown  that  the  pressure  correction  equation  can  be  derived  either  by  using  the  geometric 
conservation  equation  or  the  overall  mass  conservation  equation.  Depending  on  which 
equation  is  used,  the  segregated  pressure-based  multi-fluid  flow  algorithms  were  classified 
respectively  as  either  the  Geometric  Conservation  Based  family  of  Algorithms  (GCBA)  or 
the  Mass  Conservation  Based  family  of  Algorithms  (MCBA).  Moukalled  et  al.  [43-46] 
implemented  and  tested  the  MCBA  family  and  proved  its  capability  to  predict  multi-fluid 
flow  at  all  speeds.  On  the  other  hand,  the  GCBA  family  has  not  yet  been  implemented  nor 
tested. 

The  objective  of  the  present  work  is  to  implement  and  test  the  GCBA  family  within  a 
structured  finite-volume  framework  with  the  convection  terms  along  the  control  volume  faces 
evaluated  using  a  High  Resolution  (HR)  scheme  applied  within  the  context  of  the  Normalized 
Variable  and  Space  Formulation  methodology  (NVSF)  [15].  To  reduce  the  overall 
computational  cost,  the  convergence  rate  is  accelerated  through  the  use  of  a  non-linear  full 
multi-grid  method.  The  discretization  scheme  is  second-order  accurate  in  space  and  first 
order  accurate  in  time. 

In  what  follows,  the  governing  equations  are  first  introduced,  followed  by  a  brief  description 
of  the  discretization  procedure.  Then  the  GCBA  algorithms  are  presented,  their  capabilities  to 
predict  multi-fluid  flow  phenomena  at  all  speeds  demonstrated,  and  their  performance 
characteristics  (in  terms  of  convergence  history  and  speed)  assessed.  For  that  purpose,  a  total 


of  twelve  laminar  and  turbulent  incompressible  and  compressible  problems  encompassing 


me  Kjeumeiric  l- unservauun  ousea  /ugurunms  jor  iviuui-r  tuia  now  ai  /ui  dpeeas 


/ 


dilute  and  dense  gas-solid,  and  bubbly  flows  in  the  subsonic,  transonic  and  supersonic 
regimes  are  solved.  In  addition,  the  performance  of  these  algorithms  is  evaluated  using  (i)  a 
single  grid  approach  (SG),  (ii)  a  prolongation  only  approach  (PG)  whereby  the  solution 
moves  in  one  direction  starting  on  the  coarse  grid  and  ending  on  the  finest  grid  with  the 
solution  obtained  on  level  n  used  as  initial  guess  for  the  solution  on  level  (n+1),  and  (iii) 
finally  a  Full  Multi-Grid  (FMG)  approach  with  a  W  cycle. 

The  Governing  Equations 


In  multi-phase  flow  the  various  fluids/phases  coexist  with  different  concentrations  at  different 
locations  in  the  flow  domain  and  move  with  unequal  velocities.  Thus,  the  equations 
governing  multi-phase  flows  are  the  conservation  laws  of  mass,  momentum,  and  energy  for 
each  individual  fluid.  For  turbulent  multi-phase  flow  situations,  an  additional  set  of  equations 
may  be  needed  depending  on  the  turbulence  model  used.  These  equations  should  be 
supplemented  by  a  set  of  auxiliary  relations.  The  various  conservation  equations  needed  are: 
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at 
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where  the  meanings  of  the  various  terms  are  as  given  in  the  nomenclature. 

The  effect  of  turbulence  on  interfacial  mass,  momentum,  and  energy  transfer  is  difficult  to 
model  and  is  still  an  active  area  of  research.  Similar  to  single-fluid  flow,  researchers  have 
advertised  several  flow-dependent  models  to  describe  turbulence.  These  models  vary  in 


complexity  from  simple  algebraic  [47]  models  to  state-of-the-art  Reynolds-stress  [48] 
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models.  However,  the  one  used  here  is  the  two-equation  k-£  model  [49]  described  next.  The 
phasic  conservation  equations  governing  the  turbulence  kinetic  energy  (k)  and  turbulence 
dissipation  rate  (£)  for  the  kth  fluid  are  given  by: 

1  P  k  4v.[yw)  =  v,  r'^Vk®  +r|1|p<ll(G<k)-ell,)+I«  (4) 


3(r|k>0,k¥k,'>  (  i,<k>  1 

V  K  /  I  y7  /'r(k'>p(k)u(k)^(k)\=  ^7  r(k)iL_y£(k)  + 

at  V  H 


r(k)P(k)T^(cieG(k)-c2£e(k))  +  I^k) 


where  Ifkk)  and  I |.k)  represent  the  interfacial  turbulence  terms.  The  turbulent  viscosity  is 


calculated  as: 


(k) 


For  two-phase  flow,  several  extensions  of  the  k-£  model  that  are  based  on  calculating  the 
turbulent  viscosity  by  solving  the  k  and  £  equations  for  the  carrier  or  continuous  phase  only 
have  been  proposed  in  the  literature  [50-55].  In  a  recent  article,  Cokljat  and  Ivanov  [49] 
presented  a  phase  coupled  k-£  turbulence  model,  intended  for  the  cases  where  a  non-dilute 
secondary  phase  is  present,  in  which  the  k-£  transport  equations  for  all  phases  are  solved. 
Since  the  method  is  still  not  well  developed,  the  first  approach  in  which  only  the  k  and  £ 
equations  for  the  carrier  phase  are  solved  is  adopted  in  this  work.  Details  regarding  the 
specific  model  used  will  be  presented  as  needed. 

If  a  typical  representative  variable  associated  with  phase  (k)  is  denoted  by  (f>(k),  equations  (1)- 
(5)  can  be  presented  via  the  following  general  phasic  equation: 


afr<kykyk^ 

l  P  ^  yv.(r(klplk)u|kVk))=V.(r(k)r,k)Vf))+r(k,Q(k) 

ot  v  7  v  7 

where  the  expression  for  r(k)  and  Q(k)  can  be  deduced  from  the  parent  equations. 


(7) 
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The  above  set  of  differential  equations  has  to  be  solved  in  conjunction  with  constraints  on 
certain  variables  represented  by  algebraic  relations.  These  auxiliary  relations  include  the 
equations  of  state,  the  geometric  conservation  equation,  and  the  interfacial  mass,  momentum, 
energy,  and  turbulence  energy  transfers. 

Physically,  the  geometric  conservation  equation  is  a  statement  indicating  that  the  sum  of 
volumes  occupied  by  the  different  fluids  within  a  cell  is  equal  to  the  volume  of  the  cell 
containing  the  fluids. 

£r(k)  =  l  (8) 

k 

For  a  compressible  multi-phase  flow,  auxiliary  equations  of  state  relating  density  to  pressure 
and  temperature  are  needed.  For  the  kth  phase,  such  an  equation  can  be  written  as: 

pW=pw(p,rw)  (9) 

Several  models  have  been  developed  for  computing  the  interfacial  mass,  momentum,  energy, 
and  turbulence  energy  transfers  terms.  The  closures  used  in  this  work  will  be  detailed 
whenever  they  arise  while  solving  problems. 

In  order  to  present  a  complete  mathematical  problem,  thermodynamic  relations  may  be 
needed  and  initial  and  boundary  conditions  should  supplement  the  above  equations. 

Discretization  Procedure 


The  general  conservation  equation  (7)  is  integrated  over  a  finite  volume  to  yield: 


dfr(k,o(k,(b(k)'l 

JJ  ^  ^  +  V.(r(k,p(k,u<kVk)  )1Q 


(10) 


=  ||  V.(r(k)r(k)V4>,k))d^  +||  r(k)Q(k)dft 


a  a 

Where  Q.  is  the  volume  of  the  control  cell  (Fig.  1(a)).  Using  the  divergence  theorem  to 
transform  the  volume  integral  into  a  surface  integral  and  then  replacing  the  surface  integral 
by  a  summation  of  the  fluxes  over  the  sides  of  the  control  volume,  equation  (10)  is 


transformed  to: 
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a(r'1'p!l<rn)+  i(jr+jr)=r™Q“>n 

Ol  nh=R  w  n  c  t  h 


(11) 


where  J  D  and  J  are  the  diffusive  and  convective  fluxes,  respectively.  The 


nb=e,w,  n,s,t,b 

r  (k)C 


discretization  of  the  diffusion  term  is  second  order  accurate  and  follows  the  derivations 
presented  in  [35].  For  the  convective  terms,  the  High  Resolution  SMART  [13]  scheme  is 
employed,  even  for  the  calculation  of  interface  densities,  and  applied  within  the  context  of 
the  NVSF  methodology  [15].  Substituting  the  face  values  by  their  functional  relationship 
relating  to  the  node  values  of  (f),  Eq.  (11)  is  transformed  after  some  algebraic  manipulations 
into  the  following  discretized  equation: 


A'X1'  =£a“C+b 


00 

p 


(12) 


where  the  coefficients  Apk)  and  A  ^  depend  on  the  selected  scheme  and  Bpk)  is  the  source 
term  of  the  discretized  equation  .  In  compact  form,  the  above  equation  can  be  written  as 

lA-C  +  B'p1 

.t=hp  [<r]=  —  <i3) 

Ap 

The  discretization  procedure  for  the  momentum  equation  yields  an  algebraic  equation  of  the 
form: 

=Hp[u(k>J-r(k)Dl,k,Vp(p)  (14) 

On  the  other  hand,  the  phasic  mass-conservation  equation  (Eq.  (1))  can  be  either  viewed  as  a 

phasic  volume  fraction  equation: 

rpk)  =  Hp  [r(k)  ]  (15) 

or  as  a  the  following  phasic  continuity  equation  to  be  used  in  deriving  the  pressure  correction 

equation: 


fetK’)-(r<1)p«)' 

St 


Old 

— 0+Ap[r(k)p,k)u(k).s]  =  r(k)M,k) 


where  the  ?  operator  represents  the  following  operation: 


Ap[©]  =  £®f 

f=  NB(P) 


(16) 


(17) 
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Geometric  Conservation  Based  Algorithms  (GCBA) 

The  numbers  of  equations  describing  an  n-fluid  flow  situation  are:  n  momentum  equations,  n 
volume  fraction  (or  mass  conservation)  equations,  a  geometric  conservation  equation,  and  for 
the  case  of  a  compressible  flow  an  additional  n  auxiliary  pressure-density  relations. 
Moreover,  the  variables  involved  are  the  n  velocity  vectors,  the  n  volume  fractions,  the 
pressure  field,  and  for  a  compressible  flow  an  additional  n  unknown  density  fields.  It  is  clear 
that  the  n-velocity  fields  are  associated  with  the  n-momentum  equations,  i.e.  the  momentum 
equations  can  be  used  directly  to  calculate  the  velocity  fields.  The  volume  fractions  could 
arguably  be  calculated  from  the  volume  fraction  equations,  which  means  that  the  remaining 
equation  i.e.  the  geometric  conservation  equation  (the  volume  fractions  sum  to  1)  has  to  be 
used  in  deriving  the  pressure  equation,  or  equivalently  the  pressure  correction  equation.  This 
results  in  what  is  called  here  the  Geometric  Conservation  Based  Algorithm  (GCBA). 

The  sequence  of  events  in  the  Geometric  Conservation  Based  Algorithm  (GCBA)  is  as 
follows: 

•  Solve  the  individual  mass  conservation  equations  for  volume  fractions. 

•  Solve  the  momentum  equations  for  velocities. 

•  Solve  the  pressure  correction  equation. 

•  Correct  velocity,  volume  fraction,  density,  and  pressure  fields. 

•  Solve  the  individual  energy  equations. 

•  Return  to  the  first  step  and  repeat  until  convergence. 

The  GCBA  uses  the  momentum  equations  for  a  first  estimate  of  velocities.  However,  the 
volume  fractions  are  calculated  without  enforcing  the  geometric  conservation  equation. 
Hence,  the  mass  conservation  equations  of  all  fluids  are  used  to  calculate  the  volume 
fractions.  As  such,  the  pressure  correction  equation  should  be  based  on  the  geometric 
conservation  equation  and  used  to  restore  the  imbalance  of  volume  fractions.  The  errors  in 
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the  calculated  volume  fractions  are  expressed  in  terms  of  pressure  correction  (P),  which  is 
also  used  to  adjust  the  velocity  and  density  fields. 

The  Pressure  Correction  Equation 

After  solving  the  continuity  equations  for  the  volume  fraction  fields  and  the  momentum 
equations  for  the  velocity  fields,  the  next  step  is  to  correct  the  various  fields  such  that  the 
volume  fraction  fields  satisfy  the  compatibility  equation  and  the  velocity  and  pressure  fields 
satisfy  the  continuity  equations.  For  that  purpose,  a  guess-and-correct  scheme  is  adopted. 
Correction  is  obtained  by  solving  a  pressure  correction  equation  derived  from  the  geometric 
conservation  equation.  To  start  the  derivation,  it  is  noticed  that  initially  the  volume  fraction 
fields  denoted  by  r(k)* ,  do  not  satisfy  the  compatibility  equation  and  a  discrepancy  exists  i.e. 
RESGP  =l-£r'k)*  (18) 

k 

A  change  to  r(k>*  is  sought  that  would  restore  the  balance.  The  corrected  r  value,  denoted  by 
r(k)  (r(k)  =r(k)  +r(k) ),  is  such  that 

L(r<k)')=l-L(r(k,1  =  RESGP  (19) 

k  k 

Correction  to  the  volume  fraction,  r<k)  ,  will  be  associated  with  a  correction  to  the  velocity, 
density,  and  pressure  fields,  u(k)  ,  p(k)  ,  and  P'  respectively.  Thus,  the  corrected  fields  are 
given  as: 


r(k)  =  r(k)*  +  r(k)' ,  P  =  F  +  P',  u(k)  =  u(k)*  +u(k/,  p(k)  =  p(k)°  +  p(k)' 

The  discretized  form  of  the  corrected  continuity  equation  of  phase  (k)  can  be  written  as 


(20) 


St  P  (21) 

+  Ap((r<k)*  +r(k)  )(p(k)°  +  p(k>  )(u<k)*  +u(k)').s)  =  M(pk)(r](k)*  +4k)>p 
Neglecting  second  and  third  order  terms  (i.e.  rpk)ppk) ,  p(pk’u(k,,rpk,u(k)  ,and  ipk)  p')' ulk)  ),  its 


expanded  form  reduces  to: 
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i  j> 


(rrpj.ty+W) 

8t 


+Ap[(r(k)*p(k)Vk)'.S  +  r(k)*U(k)  pW'  +p(k)°U(k)*r(k)')] 


f  <k)*  (k)A_( /k) do)0 

_  jyj(k)j-(k)  ^2  —  Up  Pp  L  XP  Pp  / 

P  P  P  8t 


-  Ap[(r(k,*p(k)°U(k)*)]+  Mpk)rpk)*np 


Writing  utk)  as  a  function  of  P',  similar  to  what  is  usually  done  in  a  SIMPLE-like  algorithm, 
the  correction  momentum  equations  become 


u(k)-  =  HP[u(k)  ]  -  r(k)*D(1°VP'  - r(k)/D(k)VP°  - r(k)D(k)VP/ 

Substituting  Eq.  (23)  into  Eq.  (22),  rearranging,  and  discretizing  one  gets 

r«>--Hp|r<k>]= 

Ab-o  r  |'HP[ulkl]-r'k)'D,k'VP''| 

^-^p'k,'+Ap  r(k)*p<k)°  .S  +  r(krU(k,*p(k) 

8t  -r(k)D(k)VP/ 


-np+Ap|r(B'pll)'Ull,')]-M<Br“l'£2p 


where  Rpk>  =  1/Apk) . 


Neglecting  the  correction  to  neighboring  cells,  equation  (24)  reduces  to: 


r(k)*0 

p  £2pp(pk)  +  Ap  r(k)*p(k,° 


r(k)'_  p(k) 

rp  -  _kp 


HP[u(k)  ]  -  r(k)*D(k)VP/ 


.r(k)'D(k)VP' 


S  +  r(k)*U(k)  p<k,' 


(wtthiv?! 


-np+Ap[(r(k>-plkl‘Ulk>*)]-M<k>r«k>‘n[ 


Substituting  this  equation  into  the  geometric  conservation  equation  and  replacing  density 
correction  in  terms  of  pressure  correction  (i.e.  p(kl  =  Cj,k)P/ ),  the  pressure  correction  equation 
is  obtained  as 


r(k)*0  r(k) 
XP  “p'-p 


Pp  +  Ap[r(k)*U<k)*C‘k,P/]  + 


-RPk)  Ap[-(k)‘p(k)0(HP[u(k),]-r(k)tD(k)VP/-r(k),D(k)VP/)ls] 

k 

(r(k)»  (k)A_(  (k)  (k))01d 

+  l  p  Pp  ’  [  p  Pp  ’  Q.p  +  Ap  [(r(k)*p(k)°U(k)* )] 


=  RESG, 
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li- 


If  the  HP  [u  ]  term  in  the  above  equation  is  retained,  there  will  result  a  pressure  correction 
equation  relating  the  pressure  correction  value  at  a  point  to  all  values  in  the  domain.  To 
facilitate  implementation  and  reduce  cost,  simplifying  assumptions  related  to  this  term  have 
been  introduced.  Depending  on  these  assumptions,  different  algorithms  are  obtained.  A 
summary  of  the  various  GCBA  algorithms  (GCBA-SIMPLE,  GCBA-SIMPLEC,  GCBA- 
PISO,...)  used  in  this  work  is  given  next. 
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i 


ER 


frw* o  r(k> 

Ip  iipV-p 


St 


P;  +  Ap  [r(k)*U(k,*Cpk)P/]- Ap[r<k,*p(k)°(r(k)*D(k,VP/).s] 


E 


R 


(k) 


OmlMsAp?’)' 

St 


Did 

— n„+Ap[(, 


rltrp . U“’')] 


M 


+ 


Ap  [r(k)*p(k)°  (hP[u  (kr  ]  -  r(k)D(k)VP/).s] 


-  RESG , 


Approximation: 


Neglect:  HP[uw  ],  r(kyDwVP' 


u 


(k)'  --r(k)*D(k)VDP/ 


» p  A  p  M-r  p  "  pi 

Approximate  Equation: 


(34) 


(35) 


R 


(k) 


r  r‘k)*QpC(k) 


St 


—  Pp  +  Ap  [r(k)*U(k)*Cpk,P/]- Ap[r(k)*p<k)°(r(k)*D(k,VP,).s] 


-L  < 

k 


R 


(k) 


(r^pp-^p  r) 

St 


Did 

— ^p+Ap[(r(k)*pw°U(k)‘)] 


-RESG, 


(36) 


A  Global  GCBA-SIMPLE  Iteration 


•  Solve  implicitly  for  the  volume  fraction  fields. 

•  Solve  implicitly  for  u(k),  using  the  old  pressure ,  density,  and  volume  fraction  fields. 

•  Calculate  the  D(k)  fields. 

•  Solve  the  pressure  correction  equation. 

•  Correct  u(k),  P,  /k)  and  p(k|. 

•  Solve  implicitly  the  energy  equations  and  update  the  density  fields. 

•  Return  to  the  first  step  and  iterate  until  convergence. _ 


The  GCB A  following  SIMPLEC  (GCBA-SIMPLEC):  Symbolic  Form 
Predictor: 

4k)*  =HP[r,k)‘]  (37) 

u'k)*  =HPp[u(k)*]-4k,*D;,k)VpP0  (38) 

Corrector: 


u 


(k)* 


=  u(k)*  +u(k)',P*  =P°  +P' 


,(k)* 


p(k)°  +  p(k)  ,r(k)* 


.  r(k)*  +r(k) 


<k)**  =HPP[u,k^] 


(k)* 


„(k)* 

1P 


■  r±~'  DpkVpP’  =HPp[u1^  +u 


(k)» 


J 

(k)']-(: 


rpk)*  + 


(39) 


rpk)  )Dpk,Vp(p°  +  P')  (40) 


u(kV  =HPp[u(k)']-r<k)*D'k)VpP,-r'k,'D<k)VpP° -r'k)'D(pk)VpP' 


(41) 


Subtracting  HPp  [l]upk)  from  both  sides,  one  gets 
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iu 


Upk)  -HPP[l]upk>  =  HPP[u'k|-]-HPp[l]u«-  -r<k>-Df>VpP' 


-rpK,DpK,VpP°  —  fpk)  Dp  )VpP/ 


l-HPP[l]  |upk)  =  HPP[u(k)/  -  u  pk)<  ]  -  r  pk  }*D  pk)  V  P  P'  -  rpk  pk)  V  P  P  °  -if)D?)VPP' 


HPp[u(k)  —  Upky  ] 


l-HPp[l]  1  -  HPp  [l] 


l-HPP[l] 


V  P°  -r( 

v  pX  Ip 


1  -  HP  p  [l] 


V.P' 


M  HPp[u,k>  -Up*] 


l-HPp[l] 


-  rpk)*Dpk)VpPr  -  rpk)"Dpk)VpP°  -  if >  D<,k)VpP' 


p(k)'  =  C(k)P' 


(j^)* 

rpk)"  =  -Rpk)  ^ - -Ppk)  +  Ap  |r(k)*p(k)ou<k)  ).S  +  r(k)*u(k)*p(k)] 

St 


Condition: 

£{r«'}=RESG, 


r(k)*^  C(k) 

Ip 


P;  +  Ap  [r(k)*U(k)*Cpk)P/]+ 


.-•E  -Rf1  Ap  r(k)*p(k)° 


HP[u(k)  Upk)  ]  „<k)*fi(k)Y7r)/  „(krft(k)V7r)/ 


l-HPP[l] 


-rw  D  VP  -r  D  VP  S  r  =RESG 


(4k)tp^°)-fak)p‘k))C 

8t 


-Qp  +  Ap[(r(k)*p(k)°U(k)*)] 


I'  (  !C  )  ^  ^  ^ 

R(k,  fp - P^p_P'  +  Ap [r (k)*U(k)*Cpk)P/] - Ap [r (k)*p(k)° (r (k)*D<k) VP'jts] 

St 


Ik"  & 


■£2,+A,[(r,k’-p‘«'U>k|-)]- 


Af  [i',k’’p'k>'(HP[u,k>'  -  Up"  ]  -r‘krD<k,VP').S 


RESG, 


me  Kjeumeiric  l- unservauon  nasea  /ugorunms  jor  iviuui-r  tuia  now  til  /ui  jpeeas 


Approximation: 

Neglect:  HP[uw'-<r],  r<k)D(k)VP' 

=>  Upk)"  = -rpk)*Dpk)VpP/ 

Approximate  Equation: 


R 


(k) 


^r(k)*D 

Ip  ifip'-'p 


8t 


P;  +  Ap  [r  (k)*u(k)*Cpk)P/]  -  A  p  [r  (k)*p(k)°  (r  (k)*D<k)  VP/)lS_ 


~L  < 

k 


R 


(k) 


(rpk)*Ppk>° )~  (rpk)Pp° )' 

8t 


pid 

—  Qp+ApIr^p^U^)] 


-  RESG , 


A  Global  GCBA-SIMPLEC  Iteration 


(49) 


(50) 


•  Solve  implicitly  for  the  volume  fraction  fields. 

•  Solve  implicitly  for  u(k),  using  the  old  pressure ,  density,  and  volume  fraction  fields. 

•  Calculate  the  D(k)  fields. 

•  Solve  the  pressure  correction  equation. 

•  Correct  u(k),  P,  /k)  and  p(k|. 

•  Solve  implicitly  the  energy  equations  and  update  the  density  fields. 

•  Return  to  the  first  step  and  iterate  until  convergence. _ 


The  GCB A  following  PRIME  ( GCBA-PRIME):  Symbolic  Form 
Predictor: 

rpk)*  =  Hp[r(k)°]  (51) 

<k)*  =HPP[u(k)o]-r‘k),D'k,VpP°  (52) 

Corrector: 


(u'k>  ,P’,p<k>  ,rtk) 


u 


(k)* 


:U(k)*  +U(k)',P*  =P°  +  P', 


p(k)*  _  p(k)°  _|_p(k)'  J.(k)**  _  r(k)» 


(53) 


u 


(k)* 


:  HPp[u(k)**]  — rpk)**Dpk)VpP*  =  HPp[u(k)*  +  u(k)]  -  (fpk)*  +  r£k)<)Dpk)Vp(p°  +  P')  (54) 


<r  =HPp[u(k)*-u(k)o]  +  HPp[u,k,']-r1(k)*D[)k'VpP,-if,'D[)k)VpP0  -r^D^VpP' 

P(kr  =  C(k)P'  (55) 


„(k)'  _  _p(k) 
-  P  -tv  P 


(  *-(k)* 


8t 


^kr 


1  +  Ap[(r(k)*p(k)ou<k)  ).S  +  r(k)*U(k)*p(k)  ] 


Condition: 

£{r«}=RESGp 


k 


(56) 
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-R 


r(  C(  r  i 

P  P  p  -P;  +  Ap[r(k)*U(k)*C‘k,P/] 


8t 


+ 


IR 


(k) 


+ 

V 

(r-W 


r(k)yk)° 


HP[u (k)*  -u(k)o]  +  HP[u(k)  ] 

-r(k)*D(k)VP'  -  r(k)'D(k)VP' 

(c’prHrtiT 


St 


■np+Ap[(r<k»VB'U»)] 


RESGt 


R 


ip  iip'-p 

St 


(k) 


P;  +  Ap  [r(k)*U(k,*Cpk)P/]  -  Ap[r(k,*p(k)°  (r(k)*D(k,VP').s] 


(k) 


frrprH^K’J 

St 


tOld 

—  ^p+Ap[(r(kVk,°U(k,*)]+ 


Ap  [r(k)*p(k)°  (HP[u(k)*  -  u(k)°  ]  +  HP[u(k)  ]  -  r(k,D(k)  VP').s] 
Approximation: 

Neglect:  HP[u(k)*  -u(k)o],  HP[u(k)'],r(k)'D(k)VP' 

=>  u'k)'  =-r'k)tD'k,VpP/ 

Approximate  Equation: 

vl!"  LIS' 


RESG , 


R 


(k) 


St 


^P;+Ap[r(k)SU(k)*C;k,P,]-Ap[r(k,*p(k)°(r(k),D,k,VP/).s] 


R 


(k) 


St 


-£lp+Ap[(r|1>-pll>-U<1,‘)] 


M 


RESG, 


(57) 


(58) 


(59) 


(60) 


A  Global  GCBA-PRIME  Iteration 


•  Solve  explicitly  for  the  volume  fraction  fields. 

•  Solve  explicitly  for  u(k),  using  the  old  pressure  and  density  fields. 

•  Calculate  the  D(k)  fields. 

•  Solve  the  pressure  correction  equation. 

•  Correct  u(k),  P,  /k)  and  p(k|. 

•  Solve  implicitly  the  energy  equations  and  update  the  density  fields. 

•  Return  to  the  first  step  and  iterate  until  convergence. _ 


The  GCB A  following  SIMPLEST  (GCBA-SIMPLEST):  Symbolic  Form 
Predictor: 

=  H°  [r(k)*]  +  Hp  [r(k)o] 

u<k)*  =  HPpD  [u ,k)*  ]  +  HPpc  [u(k)°  J-  rpk)*Dpk)VpP° 


(61) 

(62) 
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i  y 


Corrector: 


(  (k)'  r>'  _  (k)'  ( 

lu  ,P,p  ,r 


u(k)‘*  =  u(k)*+u(k)',P*  =F+P', 
p(k)*  =  p(k)°  +p(k)',r(k)“  =  r(k)*  +r(k,/ 


=  HP°[u(k)**]+HPp  [u(k)**]-rpk)**Dpk)VpP* 

=  HP,?[u(k)*  +u(k)  ]+HPp  [u(k)*  +u(k)< ]-(rpk)*  +rpk)  )Dpk)Vp(p°  +  P') 
=  HPpD  [u(k)*]+  HPpD  [u,k/]+  HPp  [u(k)* ]  +  HPpc  [u(k/] 

-rpk)*Dpk)VpP°  -r^^Dp^VpP"  -ipk>Dpk,VpP°  -r^D^VpP' 


u<k)  =  HPp[u(k)  ]  +  HPp  [u(k)*  -  u(k)o]  -  rpk)*Dpk)VpP/-rpk)Dpk)VpP°  -  r"  D[,k,VpP 


p(k)'  =C(pk)P' 

(  Ek>*o 

«(k)  _  n  (k)  rp  ^P 


)'kV  +  Ap  [(r(k)*p  (k)ou(kr  \S  +  r«*U(k)*p(k)'  ] 


Condition: 

£{r«}=RESG, 


r(k)*o  r(k)  r 

"  "  p  Pp  +Ap[r(k,‘U(kr 

8t 


C'k,P']- 


,,,  ,ky/ HP[u  k)  ]  +  HPc[u  k  -u(k,0]^| 

V  -R‘k)  Ap  r  k  p  k)  ,S 

V  -  r(k)*D(k)VP/-  r(k)D(k)VP' 


(  (k)*0(kr)_(  (k)0(k))old 

+  lp  Pp  ’  l  p  Pp  ’  £>p  +  Ap  [(r(k)*p(k) 

8t 


RESG , 


■  u"  ■■ )] 


Er 


r(k) 

,k)  p  p;  +  Ap  [r  <k)*U  (k,*C  pk)Pr]  -  A  p  [r <k)*  p(k)°  (r  <k)*D (k)  VP').s] 


f  fr(k)*o(k)° )— (r(k,D<k)  )°ld  \!  m 

fe  Pp  Ubl  Pi-  l  np+Ap  (r"-pWU«»-  + 

-2^1KP  8t 

AP  [r(k)*p(k,°  (HPC  [u(k>*  -  u(k,°  ]  +  HP[u(k)  ]  -  r <k)  D(k) VP')ls] 


-RESG, 


Approximation: 


Neglect:  HPc[u(k)*  -  u(k)o],  HP[uK|  ],rwDwVP 


(k)'n  Jk)'n(k)nn' 


=>  <y  =-r'k)*D‘k)VpP' 
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Approximate  Equation: 


(70) 


A  Global  GCBA-SIMPLEST  Iteration 

•  Solve  for  r,k),  treating  HD[r(k)J  implicitly  and  Hc[rlk)J  explicitly. 

•  Solve  for  u(k>,  treating  HP  D  [u( k  1 J  implicitly  and  HP 1  [u !  1 '  J  explicitly. 

•  Calculate  the  D(k)  fields. 

•  Solve  the  pressure  correction  equation. 

•  Correct  u(k),  P,  /k)  and  p(k|. 

•  Solve  implicitly  the  energy  equations  and  update  the  density  fields. 

•  Return  to  the  first  step  and  iterate  until  convergence. _ 
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The  GCBA  following  PISO  (GCBA-PISO):  Symbolic  Form 
First  Predictor: 


.(k)« 


Hp[r(k)*] 


Up  =  HPp[u  ]-r|,  D|,  VpP" 

First  Corrector: 
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Approximation: 


Neglect:  HPfu®'], r(k)D(k)VP' 

=*  u'k)'  =-r'k)*D'k,VpP/ 

Approximate  Equation: 
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Second  Predictor: 
rpk  =  Hp  [r(k)  ] 

u£)m  =  HP;*[u(k)"]  -r<k)t“D(pk)"VpP* 

Second  Corrector: 
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Condition: 
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Approximation: 

Neglect:  HP  |uik'"'  - u(k)“  +  u(k)'],  r(k)'D(k)**VP// 

=>  u  pk =  -ipk)oDpk,**VpP" 

Approximate  Equation: 
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A  Global  GCBA-PISO  Iteration 
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•  Solve  implicitly  for  /k\ 

•  Solve  implicitly  for  u  using  the  old  pressure  and  density  fields. 

•  Calculate  the  D(k)  fields. 

•  Solve  the  pressure  correction  equation. 

•  Correct  ua),  P,  and  pa  \ 

•  Solve  implicitly  the  energy  equation  and  update  the  density  fields. 

•  Solve  explicitly  for  /k). 

•  Solve  the  momentum  equations  explicitly  and  calculate  the  D(k)  fields. 

•  Solve  the  pressure  correction  equation. 

•  Correct  ua),  P ,  and  pa  \ 

•  Return  to  step  one  and  iterate  until  convergence _ 


The  GCB A  following  SIMPLEX  (GCBA-SIMPLEX):  Symbolic  Form 
Predictor: 

=  Hp[r(k)*] 
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Approximation: 


Neglect:  rpk)  Dpk)VpP',  and  let 
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(99) 
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=>  -  rpk  }*D  pk)  sx  V  pP '  =  HPp[-r(k)*D(k)SXVP/]  -  rpk)*Dpk)VpP/ 

Assume  that  the  pressure  difference  local  to  the  velocity  is  representative  of  all 
differences  i.e.  HPP [-r(k)*D(k)SXVP/]  =  -(VpP/)HPp[r(k)*D(k)SX  ] ,  thus: 
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Approximate  Equation: 
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A  Global  GCBA-SIMPLEX  Iteration 


•  Solve  implicitly  for  hk>. 

•  Solve  Implicitly  for  u(k),  using  the  old  pressure  and  density  fields. 

•  Calculate  the  D(k)  fields. 

•  Solve  implicitly  for  the  DlklSX  fields. 

•  Solve  the  pressure  correction  equation  using  these  Dlk,sx  fields. 

•  Correct  u(k),  P,  /k),  and  p(k). 

•  Solve  implicitly  the  energy  equations  and  update  the  density  fields. 

•  Return  to  the  first  step  and  iterate  until  convergence. _ 


The  GCB A  following  SIMPLEM  (GCBA-SIMPLEM):  Symbolic  Form 
First  Predictor: 

r'k)‘=HP[r(k,‘] 

Calculate  the  coefficients  of  the  momentum  equations. 

First  Corrector: 
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Approximation: 

Neglect:  HP[u(k)'],  r(ky D(k)VP' 

=>  u<,ky  =-r<k)*D<,k)VpP/ 

Approximate  Equation: 
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Second  Predictor: 


(HI) 


(112) 


(113) 


u|,k)"  =  HPp  [u (k)** ]  -rpk)“Dpk)*VpP* 

Second  Corrector: 


(114) 


No  corrector  stage. 

A  Global  MCBA-SIMPLEM  Iteration 
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•  Solve  implicitly  for  /k\ 

•  Calculate  the  Dlk)  fields  based  on  values  from  the  previous  iteration. 

•  Solve  the  pressure  correction  equation. 

•  Correct  u(k),  P,  /k),  and  p(k). 

•  Calculate  new  HP<k)  and  D(k)  fields. 

(k) 

•  Solve  implicitly  for  u  using  the  new  fields. 

•  Solve  implicitly  the  energy  equations  and  update  the  density  fields. 

•  Return  to  the  first  step  and  iterate  until  convergence. _ 


The  Expanded  Form  of  the  Pressure-Correction  Equation 


If  r(k),  U,k)  and  p(k)  denote  values  from  the  previous  iteration  or  from  a  previous  corrector  step, 
the  pressure  correction  equation,  applicable  to  all  algorithms,  becomes 
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The  discretization  of  the  above  equation  yields 
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Following  the  calculation  of  the  pressure  correction  field,  Upk)  ,  p(pk) ,  and  rpk)  are  obtained 


using  the  following  equations 
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The  Multi-Grid  Strategy 


(118) 


Similar  to  other  iterative  methods,  the  rate  of  convergence  of  the  solution  method  described 
above  does  not  scale  linearly  with  the  grid  size,  rather  the  convergence  rate  decreases  more 
drastically  as  the  number  of  grid  points  increases.  This  behavior  is  attributed  to  the  speed  at 
which  the  iterative  solver  transports  the  boundary  information  across  the  domain  (e.g.  with 
SOR  one  grid  point  per  iteration).  Since  information  has  to  travel  back  and  forth  several 
times  to  achieve  convergence,  acceleration  of  the  outer  iterations  through  the  use  of  multi¬ 
grid  methods  is  essential  when  solving  over  large  grids.  The  idea  underlying  the  multi-grid 
strategy  is  to  use  progressively  coarser  grids  to  accelerate  the  convergence  rate.  In 
mathematical  terms  the  low-frequency  error  components  in  the  finest  grid  appear  on  coarser 
grids  as  high-frequency  Fourier  mode  that  can  be  resolved  efficiently  by  iterative  relaxation 
solvers.  In  the  present  work,  this  strategy  is  adopted  to  accelerate  convergence  and  thereby 
reduce  the  overall  computational  cost.  The  method  used  is  the  FMG-FAS  method  [56].  For  a 
review  of  Multi- grid  methods  the  reader  is  referred,  among  others,  to  [56,57],  therefore  it  is 
sufficient  here  to  give  a  general  description  of  the  method  used. 

The  multi-grid  algorithm  adopted  in  this  work  can  be  summarized  as  follows.  Starting  with 
the  fine  mesh,  the  coarser  grid  cells  are  generated  through  agglomeration  of  four  finer  grid 
cells,  two  in  each  direction.  On  the  other  hand,  if  a  finer  grid  is  required,  subdividing  the 
coarser  grid  control  volume  into  four  control  volumes,  again  two  in  each  direction,  generates 
its  control  volumes.  With  the  FMG  cycle,  the  algorithm  starts  at  the  coarsest  level,  where  the 
solution  is  first  computed;  this  solution  is  interpolated  onto  the  next  finer  mesh,  where  it  is 


me  Kjeumeiric  l- unservauun  nasea  /ugorunms  jor  iviuiu-r  tuia  now  til  /ui  dpeeas 


used  as  initial  guess.  This  stage  is  called  the  prolongation  stage  (see  Fig.  1(b)).  Then 
iterations  are  performed  on  the  fine  mesh  and  the  solution  is  transferred  back  to  the  coarser 
mesh  by  applying  a  restriction  operator.  In  order  to  obtain  the  same  approximation  on  each 
level,  a  forcing  term  is  added  to  the  discrete  conservation  equations  on  the  coarser  grid.  This 
term  represents  the  truncation  error  on  the  coarse  grid  with  respect  to  the  fine  grid.  After 
performing  a  number  of  iterations  on  the  coarse  mesh,  the  solution  is  transferred  back  to  the 
finer  mesh  in  the  form  of  a  correction  and  a  number  of  iterations  are  performed  on  the  finer 
grid  to  smooth  the  fields.  This  process  is  continued  until  a  converged  solution  on  the  fine 
mesh  is  obtained  (see  Fig.  1(c)).  Then  the  solution  is  extrapolated  to  correct  the  finer  mesh 
fields,  followed  by  a  number  of  smoother  iterations  on  the  finer  mesh  and  the  process 
repeated  until  convergence  is  reached  on  the  desired  finest  mesh.  This  strategy  has  been 
applied  to  both  incompressible  and  compressible  supersonic  multi-fluid  flows  and  good 
savings  have  been  realized  as  will  be  shown  in  the  results  section. 

In  the  restriction  step  the  coarse  grid  variables  are  computed  from  the  fine  grid  values  as: 

$  C  =  X  E  4  F.  +  V  4>  Fi  '  d  FiC  )  (119) 

while  in  the  prolongation  step  the  fine  grid  corrections  are  computed  from  the  coarse  grid 
values  as 

<KFl  =  <I>C  +  V(KC  •  dCFi  (120) 

where  d  is  the  position  vector  connecting  points  C  and  F,  and  4k  given  by 

4>c=4>c-ic  (121) 

The  special  character  of  the  volume  fraction  and  k-8  equations  necessitates  modification  to 

the  prolongation  procedure  as  described  next. 

While  extrapolating  the  volume  fraction  field  from  the  coarse  to  the  fine  grid,  the 
prolongation  operator  may  yield  negative  volume  fraction  values  or  values  that  are  greater 
than  one.  Such  unphysical  values  are  detrimental  to  the  overall  convergence  rate  and  may 
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z, y 


cause  divergence.  To  circumvent  this  problem,  a  simple  yet  very  effective  treatment  is 
adopted:  Once  the  r-values  are  extrapolated,  a  check  is  performed  to  make  sure  they  are 
within  bounds.  If  any  of  the  r-values  is  found  to  be  unbounded,  the  r-phasic  volume  fraction 
equation  is  solved  starting  with  the  interpolated  values  until  all  of  the  r-values  are  within  the 
set  bounds.  Typically  less  than  10  iterations  are  needed.  This  treatment  has  been  found  to  be 
very  effective  and  to  preserve  the  convergence  acceleration  rate.  The  practices  of  solving  the 
volume  fraction  equations  only  on  the  fine  grid  or  forcing  the  extrapolated  unbounded  values 
to  be  within  the  set  bounds  or  discarding  corrections  that  result  in  unbounded  values  [57] 
proved  to  be  ineffective  and  slowed  the  convergence  rate  considerably. 

For  the  k-e  turbulence  model,  the  treatment  suggested  by  Cornelius  et  al.  [58]  is  adopted. 
This  approach  is  based  on  the  observation  that  the  application  of  wall  functions  to  the  coarse 
grids  would  lead  to  unphysical  values  because  of  the  relatively  large  distance  between  the 
wall  and  the  boundary  cell  center.  Thus,  at  wall  boundaries  the  restricted  fine  grid  values  of  k 
and  £  are  held  constant,  and  hence  no  corrections  are  calculated.  In  order  to  satisfy  the 
realizability  constraint,  the  restricted  turbulence  properties  and  prolongated  correction  values 
are  modified  accordingly. 

In  addition  to  the  FMG  strategy,  the  PG  approach  is  also  tested.  This  approach  differs  from 
the  FMG  method  in  that  the  solution  moves  in  one  direction  from  the  coarse  to  the  fine  grids 
with  the  initial  guess  on  level  n+1  obtained  by  interpolation  from  the  converged  solution  on 
level  n  (Fig.  1(b)).  As  such,  the  acceleration  over  the  SG  method  obtained  with  this  approach 
is  an  indication  of  the  effect  of  initial  guess  on  convergence. 

Results  and  Discussion 

The  performance  of  the  various  multi-fluid  Geometric  Conservation  Based  Algorithms  is 
assessed  in  this  section  by  presenting  solutions  to  several  one  and  two-dimensional  two-phase 
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flow  problems.  Results  are  presented  in  terms  of  the  CPU-time  needed  to  converge  the 
solution  to  a  set  level  and  of  the  convergence  history.  Moreover,  solutions  are  obtained  for  a 
number  of  grids  in  order  to  assess  the  performance  of  the  various  algorithms  with  increasing 
grid  density.  For  the  two-dimensional  problems,  in  addition  to  the  CPU-time  needed  to  solve 
a  problem  using  the  single  grid  method  (SG),  the  CPU-time  needed  using  different  solution 
strategies  is  also  displayed.  These  include  the  prolongation  scheme  (PG)  in  which  the 
solution  from  the  next  coarser  grid  is  used  to  provide  the  initial  field,  and  the  full  multi-grid 
method  (FMG).  Results  are  compared  against  available  experimental  data  and/or 
numerical/theoretical  values.  The  residual  of  a  variable  (f>  at  the  end  of  an  outer  iteration  is 
defined  as: 

RES/’^IXC1-  Ea^V-B®  (122) 

c.v  all  p  neighbours 

For  global  mass  conservation,  the  imbalance  in  mass  is  defined  as: 

Resc  =  ££ 

k  c.v. 

All  residuals  are  normalized  by  their  respective  inlet  fluxes.  Computations  are  terminated 
when  the  maximum  normalized  residual  of  all  variables,  drops  below  a  very  small  number  £s. 
For  a  given  problem,  the  same  value  of  £s  is  used  with  all  algorithms.  In  general,  it  is  found 
that  requiring  the  overall  mass  residuals  to  be  satisfied  to  within  £s  is  a  very  stringent 
requirement  and  the  last  to  be  fulfilled.  This  is  why  these  residuals  are  the  ones  presented 
here  and  used  to  compare  the  performance  of  the  various  algorithms.  In  all  problems,  the  first 
phase  represents  the  continuous  phase  (denoted  by  a  superscript  (c)),  which  must  be  fluid, 
and  the  second  phase  is  the  disperse  phase  (denoted  by  a  superscript  (d)),  which  may  be  solid 
or  fluid.  Unless  otherwise  specified  the  HR  SMART  scheme  is  used  in  all  computations 
reported  in  this  study.  For  a  given  problem,  all  results  are  generated  starting  from  the  same 
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initial  guess.  Moreover,  it  should  be  stated  that  in  iterative  techniques,  different  initial 
guesses  might  require  different  computational  efforts. 

One-dimensional  two-phase  validation  problems 

Due  to  the  large  number  of  parameters  affecting  the  performance  of  the  various  multi-phase 
Geometric  Conservation  Based  Algorithms  and  to  allow  a  thorough  testing  of  these 
algorithms,  eight  one-dimensional  two-phase  flow  problems  are  considered.  These  problems 
can  be  broadly  classified  as:  (i)  horizontal  particle  transport,  and  (ii)  vertical  particle 
transport.  Results  are  presented  in  terms  of  the  convergence  history  and  the  CPU-time  needed 
to  converge  the  solution  to  a  set  level.  Predictions  are  compared  against  available 
numerical/theoretical  values. 

Despite  its  geometric  simplicity,  the  one-dimensional  particle  transport  problem  can  represent 
a  wide  range  of  physical  conditions.  The  effects  of  grid  refinement  on  accuracy  and 
convergence  are  studied  by  solving  the  problems  on  four  grid  systems  of  sizes  20,  40,  80,  and 

Q 

160  control  volumes  with  £s  assigned  the  value  of  10"  . 

Many  runs  were  performed  so  as  to  set  the  control  parameters  of  each  algorithm  near 
optimum  values.  To  allow  a  comparative  assessment  of  performance,  the  CPU  times  are 
reported  in  the  form  of  charts.  Moreover,  all  CPU  times  are  normalized  by  the  time  needed 
by  GCBA-SIMPLE  to  reach  the  set  residuals  on  the  coarsest  grid. 

Horizontal  particle  transport 

The  physical  situation  is  depicted  in  Fig.  1(d).  Depending  on  the  set  densities,  it  represents 
either  the  steady  flow  of  solid  particles  suspended  in  a  free  stream  of  air  or  the  steady  flow  of 
air  bubbles  in  a  stream  of  water.  The  slip  between  the  phases  determines  the  drag,  which  is 
the  sole  driving  force  for  the  particle-bubble/air- water  motion  (g=0).  In  the  suspension,  the 
inter-particle/bubble  forces  are  neglected.  Diffusion  within  both  phases  is  set  to  zero  while 
the  inter-phase  drag  force  is  calculated  as: 


me  Kjeumeiric  l- unservauon  nasea  /ugorunms  jor  iviuui-r  tuia  now  til  /ui  jpeeas 


j/- 


3  C, 


!m  =  “~r<d)p(C)Vslip(u(d)  -U(c)) 


8  r. 
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r  P  XUp(u  -U  ) 


\r  (d)  (c) 

Vslip  =  U  -u 


(124) 

(125) 

(126) 

The  drag  coefficient,  Co,  is  set  to  0.44.  Since  phasic  diffusion  is  neglected,  the  GCBA- 
SIMPLEST  and  GCBA-PRIME  becomes  identical  and  reference  will  be  made  to  GCBA- 
SIMPLEST  only.  The  task  is  to  calculate  the  particle/bubble-velocity  distribution  as  a 
function  of  position.  If  the  flow  field  is  extended  far  enough  (here  computations  are 
performed  over  a  length  of  L=2m),  the  particle/bubble  and  fluid  phases  are  expected  to 
approach  an  equilibrium  velocity  given  by: 


U  =r(c)V(c)  +r(d)  V<d) 

equilibrium  inlet  inlet  inlet  inlet 


(127) 


Problem  1 :  Dilute  gas-solid  flow 


The  steady  flow  of  dilute  particles  suspended  in  a  free  stream  of  air  is  studied  first.  At  inlet, 
the  air  and  particle  velocities  are  5  m/s  and  1  m/s,  respectively.  The  physical  properties  of  the 
two  phases  are:  p(d)  /p(c)  =  2000,  r  =  1  mm,  r^t  =  10  \  Due  to  the  dilute  concentration  of 
the  particles,  the  free  stream  velocity  is  more  or  less  unaffected  by  their  presence  and  the 
equilibrium  velocity  is  nearly  equal  to  the  inlet  free  stream  velocity.  Based  on  this 
observation,  Morsi  and  Alexander  [59]  obtained  the  following  analytical  solution  for  the 
particle  velocity  u<d)  as  a  function  of  the  position  x  and  the  properties  of  the  two  phases: 


Ln[vil„t1-u«>]+— ^ 


3  pw  C, 


'  X  + 


Ln[V, 


(c) 

inlet 


■  Vtd)  1+ 

v  inlet  J  ^ 


v; 


(c) 

inlet 


(128) 


t-(c)  ..(d)  q  ~(d)  L  inlet  v  inlet  J  xr(c)  y(d) 

inlet  ~  U  8  P  rp  ’inlet  ~  Vinlet 

This  case  is  of  particular  importance  since  the  flow  situation  has  an  exact  solution.  As  shown 
in  Fig.  2(a)  the  predicted  particle  velocity  distribution  falls  on  top  of  the  analytical  solution 
given  by  Eq.  (128),  which  is  an  indication  of  the  accuracy  of  the  numerical  procedure.  The 


convergence  histories  of  the  various  GCBA  over  the  four  grid  networks  used  are  displayed  in 


Figs.  2(b)-2(h).  For  all  algorithms,  the  required  number  of  iterations  increases  as  the  grid  size 
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increases,  with  PISO  (Fig.  2(b))  requiring  the  minimum  and  SIMPLEST/PRIME  (Fig.  2(f)) 
the  maximum  number  of  iterations  on  all  grids.  The  convergence  histories  of  SIMPLE, 
SIMPLEC,  SIMPLEM,  and  SIMPLEX  (Figs.  2(c),  2(d),  2(e),  and  2(g),  respectively)  are  very 
similar  with  SIMPLEM  (Fig.  2(e))  requiring  the  lowest  number  of  iterations.  The 
convergence  paths  of  the  various  algorithms  over  a  grid  of  size  80  C.V.  are  compared  in  Fig. 
2(h)  and  the  above  observations  are  easily  inferred  from  the  figure. 

Problem  2:  Dense  gas-solid  flow 

The  only  difference  between  this  case  and  the  previous  one  is  in  the  concentration  of 
particles,  which  is  set  to  rr^t  =  10  2 .  Despite  the  low  value  of  the  inlet  disperse  phase  volume 
fraction,  the  ratio  of  disperse  phase  and  continuous  phase  mass  loadings  is  large 
r(dlp'd)/ r(c)pU)  =  20.  Thus  the  disperse  phase  carries  most  of  the  inertia  of  the  mixture.  The 
equilibrium  velocity  in  this  case,  as  obtained  from  Eq.  (127)  is  4.96  m/s  as  compared  to 
4.99996  m/s  in  the  previous  case.  Due  to  this  slight  difference  between  the  inlet  air  velocity 
and  the  final  equilibrium  velocity,  the  free  stream  velocity  may  be  assumed  to  be  nearly 
constant  and  the  variation  in  particle  velocity  can  be  obtained  again  from  Eq.  (128).  The 
predicted  air  and  particle  velocity  distributions  are  displayed  in  Fig.  3(a).  The  numerical  and 
analytical  particle  velocity  profiles  are  indistinguishable  and  fall  on  top  of  each  other. 
Moreover,  the  slight  decrease  in  the  air  velocity  can  be  easily  depicted.  The  convergence 
paths  for  all  algorithms  and  over  all  grid  systems  used  are  displayed  in  Figs.  3(b)-3(h).  In 
general,  higher  number  of  iterations  is  required  to  reach  the  desired  level  of  convergence  on  a 
given  grid  as  compared  to  the  dilute  case  due  to  the  increased  importance  of  the  inter-phase 
term.  The  general  convergence  trend  is  similar  to  that  of  the  dilute  problem  with  PISO 
requiring  the  minimum  and  SIMPLEST  the  maximum  number  of  iterations.  The  SIMPLEM 
algorithm  (Fig.  3(e))  is  seen  to  require  a  slightly  lower  number  of  iterations  on  the  finest  grid 


as  compared  to  SIMPLE  (Fig.  3(c)),  SIMPLEC  (Fig.  3(d)),  and  SIMPLEX  (Fig.  3(c)).  As 
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depicted  in  Figs.  3(f)  and  3(h),  the  performance  of  SIMPLEST/PRIME  is  poor  as  compared 
to  other  algorithms  for  the  same  reasons  stated  above. 

Problem  3:  Dilute  bubbly  flow 

For  the  same  configuration  displayed  in  Fig.  1(d),  the  continuous  phase  is  considered  to  be 
water  and  the  disperse  phase  to  be  air.  The  resulting  flow  is  denoted  in  the  literature  by 
bubbly  flow.  With  the  exception  of  p(d) /p<c)  =  10~3and  at  inlet  =0.1,  other  physical 
properties  and  inlet  conditions  are  the  same  as  those  considered  earlier.  This  is  a  strongly 
coupled  problem  and  represents  a  good  test  for  the  numerical  procedure  and  performance  of 
the  algorithms.  The  correct  physical  solution  is  that  the  bubble  and  continuous  phase 
velocities  both  reach  the  equilibrium  velocity  of  4.6  m/s  (Eq.  (127))  in  a  distance  too  small  to 
be  correctly  resolved  by  any  of  the  grid  networks  used.  Results  for  this  case  are  presented  in 
Fig.  4.  Axial  velocity  distribution  for  both  water  and  air  are  displayed  in  Fig.  4(a).  As 
expected,  both  phases  reach  the  equilibrium  velocity  of  4.6  m/s  over  a  very  short  distance 
from  the  inlet  section  and  remain  constant  afterward.  The  relative  convergence  characteristics 
of  the  various  algorithms  remain  the  same.  However,  all  algorithms  require  larger  number  of 
iterations  as  compared  to  the  dilute  gas  solid  flow  case  due  to  the  stronger  coupling  between 
the  phases.  Consistently,  the  PISO  (Fig.  4(b))  and  SIMPFEST/PRIME  (Fig.  4(f))  algorithms 
need  the  lowest  and  highest  number  of  iterations,  respectively.  As  in  the  previous  two  cases, 
the  convergence  attributes  of  SIMPFE  (Fig.  4(c)),  SIMPFEC  (Fig.  4(d)),  SIMPFEM,  and 
SIMPFEX  (Fig.  4(g))  are  very  similar  with  SIMPFEM  consistently  requiring  a  lower  number 
of  iterations.  The  large  difference  in  performance  between  SIMPFEST/PRIME  and  the 
remaining  algorithms  is  clearly  demonstrated  in  Fig.  4(h). 

Problem  4:  Dense  bubbly  flow 

The  only  difference  between  this  case  and  the  previous  one  is  in  the  concentration  of  bubbles, 
which  is  set  to  r‘n^t  =0.5 .  With  such  high  value  of  void  fraction,  bubble  coalescence  may 
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occur.  However,  this  is  not  accounted  for  here.  The  analytical  solution  is  the  same  as  in  the 
previous  case  with  the  equilibrium  velocity,  as  computed  from  Eq.  (127),  being  3  m/s.  As 
depicted  in  Fig.  5(a),  the  equilibrium  velocity  obtained  numerically  is  exact.  With  the 
exception  of  requiring  higher  number  of  iterations  to  reach  the  desired  level  of  convergence, 
the  performance  of  the  various  algorithms  (Figs.  5(b)-5(h))  vary  relatively  in  a  manner 
similar  to  what  was  previously  discussed  and  deemed  redundant  to  be  repeated. 

CPU  time:  Horizontal  particle  transport 

The  normalized  CPU  efforts  required  by  the  various  algorithms  over  all  grids  are  depicted  in 
Fig.  6.  The  charts  clearly  show  that  the  CPU  time  increases  with  increasing  grid  density.  For 
the  dilute  gas-solid  problem  (Fig.  6(a)),  it  is  hard  to  see  any  noticeable  difference  in  the  CPU 
times  for  SIMPFE,  SIMPFEM,  and  SIMPFEX.  The  SIMPFEC  and  PISO  algorithms  require 
slightly  lower  and  higher  computational  efforts,  respectively,  as  compared  to  SIMPFE.  The 
worst  performance  is  for  SIMPFEST  which  degenerates  to  PRIME  in  the  absence  of 
diffusion  and  results  in  a  fully  explicit  solution  scheme.  For  the  dense  gas-solid  flow  (Fig. 
6(b)),  the  computational  times  needed  by  SIMPFE,  SIMPFEC,  SIMPFEM,  and  SIMPFEX 
are  nearly  identical.  PISO,  however,  requires  higher  computational  effort  (50%  more  than 
SIMPFE  on  the  finest  meshes  (80  and  160  C.V.)).  The  computational  effort  needed  by 
SIMPFEST/PRIME  is  however  the  most  extensive  and  is  nearly  500%  the  one  needed  by 
SIMPFE  on  the  finest  mesh. 

The  normalized  CPU  time  of  SIMPFEST/PRIME  for  the  bubbly  flow  problems  (Figs.  6(c) 
and  6(d))  is  lower  than  in  the  previous  two  problems  due  to  a  higher  rate  of  increase  in  the 
time  needed  by  other  algorithms  (the  computational  time  of  all  algorithms  has  increased). 
The  relative  performance  of  the  various  algorithms  is  nearly  as  described  earlier  with  the  time 
required  by  of  PISO,  SIMPFE,  SIMPFEC,  and  SIMPFEX  being  on  average  the  same.  The 
SIMPFEST/PRIME  algorithm  however,  requires  nearly  50%  more  time  than  SIMPFE,  which 
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represents  a  noticeable  improvement.  The  best  performance  for  the  dense  bubbly  flow 
problem  is  for  SIMPLEM,  which  requires  about  50%  less  effort  on  the  finest  mesh  than 
SIMPLE. 

Vertical  particle  transport 

Here,  the  flow  is  in  the  vertical  direction  (Fig.  1(d)),  the  gravitational  acceleration  is  assigned 
the  constant  value  of  g=10  m/s“,  and  the  flow  field  is  extended  over  a  length  of  L=20m.  For 
this  situation,  the  velocities  of  the  two  phases  do  not  reach  an  equilibrium  value.  Rather,  the 
disperse  phase  equilibrates  to  a  finite  settling  velocity  relative  to  the  continuous  phase,  at 
which  the  gravitational  force  balances  the  drag  force  [59].  As  for  the  horizontal  transport 
problems,  the  inter-particle/bubble  forces  are  neglected.  However,  unlike  the  previous 
situation,  diffusion  in  the  continuous  phase  is  retained.  Moreover,  the  inter-phase  drag  force 
is  calculated  using  Eqs.  (124)-(126)  and  the  drag  coefficient,  Cd,  is  considered  to  be  particle 
Reynolds  number  dependent  and  calculated  as: 

CD  =  —  +  0.44,  Re  =  2rpVslip  (129) 

Re  p  p  fl c 

Since  diffusion  in  the  continuous  phase  is  not  neglected,  GCBA-SIMPLEST  and  GCBA- 
PRIME  are  expected  to  behave  differently. 

Problem  5:  Dilute  gas-solid  flow 

The  material  properties  and  boundary  conditions  considered  for  this  case  are  given  by: 

pld>/pM  =  1000,  S<c>=l(r!,  rp  =1  mm  (130) 

V“  =  lOOm/s,  V*  =10m/s,r®  =10-5  (131) 

The  large  velocity  boundary  condition  is  used  to  ensure  that  the  solid  phase  does  not  exit  the 

inlet.  The  predicted  air  and  particle  velocity  distributions  depicted  in  Fig.  7(a)  are  in  excellent 

agreement  with  similar  predictions  reported  in  [60].  As  shown  in  Figs.  7(b)-7(h),  the  mass 

residuals  tend  to  slightly  increase  at  the  beginning  of  the  iterative  process,  stagnate  over  a 

number  of  iterations  (this  number  increases  with  increasing  grid  size),  and  then  decrease 


me  Kjeumeiric  l- unservauon  nasea  /ugorunms  jor  iviuiu-r  tuia  now  til  /ui  dpeeas 


j  / 


rapidly  to  the  desired  level  of  convergence.  This  behavior  is  also  true  for  the  horizontal  case 
discussed  earlier  and  is  attributed  to  the  approximations  introduced  to  the  pressure  correction 
equation  especially  with  regard  to  neglecting  second  order  correction  terms,  which  may  be 
important  at  the  beginning  of  the  iterative  process.  Once  these  neglected  terms  become 
unimportant,  the  rate  of  convergence  increases  drastically.  Retaining  these  terms  could 
improve  the  convergence  rate  but  this  has  not  been  considered  in  this  work. 

As  depicted  in  Fig.  (7),  the  number  of  iterations  needed  to  converge  the  solution  to  the 
desired  level  is  very  close  to  that  needed  in  the  similar  horizontal  transport  case.  This  is 
equally  true  with  regard  to  the  relative  performance  of  the  various  algorithms.  The 
performance  of  SIMPLEST  (Fig.  7(f))  and  PRIME  (Fig.  7(g))  is  very  close  due  to  the  fact 
that  the  implicitness  introduced  by  the  diffusion  of  the  continuous  phase  does  not  seem  to  be 
that  important.  However,  both  require  on  the  finest  mesh  almost  430%  the  number  of 
iterations  needed  by  SIMPLE.  The  number  of  outer  iterations  needed  by  SIMPLEX  and 
SIMPLEC  is  very  close  to  that  of  SIMPLE,  while  SIMPLEM  entails  lower  number  of 
iterations.  Again,  PISO  requires  the  lowest  number  of  iterations. 

Problem  6:  Dense  gas-solid  flows 

The  material  properties  and  boundary  conditions  are  similar  to  the  previous  case  with  the 
exception  of  the  particles’  volume  fraction,  which  is  set  to  r^t  =  10~2 .  Predicted  air  and 
particle  velocity  profiles  are  displayed  in  Fig.  8(a)  while  mass  residuals  are  presented  in  Figs. 
8(b)-8(h).  Higher  number  of  iterations  is  needed  in  comparison  with  the  dilute  case  due  to  the 
higher  mass-loading  ratio.  Besides  that,  the  convergence  behavior  is  similar  to  the  previous 
cases  with  SIMPLEST  (Fig.  8(f))  and  PRIME  (Fig.  8(g))  requiring  the  highest  number  of 
iterations  and  PISO  (Fig.  8(b))  the  lowest  number  of  iterations.  The  number  of  iterations 
needed  by  SIMPLE,  SIMPLEC,  and  SIMPLEX  (compare  Figs.  8(c),  8(d),  and  8(h))  is  very 
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close.  The  SIMPLEM  algorithm  (Fig.  8(e))  needs  about  25%  less  iterations  than  SIMPLE 
(Fig.  8(c))  on  the  finest  mesh. 

Problem  7:  Dilute  bubbly  flows 

In  this  problem,  the  continuous  phase  is  water  and  the  disperse  phase  is  air.  With  the 
exception  of  p(d,/plc)  set  to  10‘3,  U(d)  =  U(c)  =1,  and  r'j’,  to  0.1  at  inlet,  other  physical 
properties  and  inlet  conditions  are  the  same  as  those  considered  earlier.  This  is  a  very 
difficult  problem  to  get  convergence  to  unless  the  proper  under-relaxation  is  used.  By  starting 
with  relatively  high  under-relaxation  factors,  the  number  of  iterations  needed  with  all 
algorithms  was  found  to  be  very  high.  In  order  to  get  feasible  solutions,  the  under-relaxation 
factors  during  the  first  20  iterations  were  set  to  0.05  and  then  increased  to  the  desired  values. 
This  was  found  to  greatly  improve  the  convergence  rate  and  to  generate  solutions  with  nearly 
the  same  computational  effort  needed  in  the  previous  cases.  In  addition,  this  treatment  has 
improved  the  performance  of  SIMPLEST  and  PRIME  dramatically  and  has  decreased  their 
required  number  of  iterations  to  values  similar  to  those  needed  by  other  algorithms  (Figs. 
9(b)-9(h)).  In  fact,  SIMPLEST  and  PRIME  are  performing  slightly  better  than  SIMPLE  for 
this  particular  problem.  Overall,  none  of  the  algorithms  shows  an  outstanding  superiority  in 
performance  over  others. 

Problem  8:  Dense  bubbly  flows 

With  the  exception  of  setting  r^  to  0.5,  the  physical  situation,  material  properties,  and 
boundary  conditions  are  the  same  as  in  the  previous  problem.  Results  for  the  problem  are 
presented  in  Fig.  10.  It  was  possible  to  get  feasible  solutions  (i.e  with  reasonable 
computational  time)  only  when  under-relaxing  by  inertia  (i.e.  through  the  use  of  false  time 
steps).  For  the  results  presented  in  Fig.  10,  a  time  step  ((it)  of  value  10‘4  s  is  used  for  the 


velocity  field  of  the  dispersed  gas  phase,  At=l  s  for  the  volume  fractions,  and  At =0.01  s  for 
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the  velocity  field  of  the  liquid  phase  for  the  grids  of  sizes  20  and  40  C.V.  and  At  =0.05  for  all 
variables  and  for  both  phases  with  the  dense  grids  (i.e.  80  and  160  C.V.).  The  predicted  liquid 
and  gas  velocity  distributions,  which  are  in  excellent  accord  with  published  data,  are  depicted 
in  Fig.  10(a).  The  trend  of  convergence  (Figs.  10(b)- 10(h))  is  different  than  what  has  been 
presented  so  far  with  all  algorithms,  except  PISO,  requiring  nearly  the  same  number  of 
iterations  and  behaving  in  almost  the  same  manner.  It  is  also  noticed  that  the  number  of 
iterations  needed  on  the  finest  mesh  is  smaller  than  the  numbers  needed  on  the  grids  of  sizes 
40  and  80  C.V.  Nevertheless,  it  was  possible  to  obtain  solutions  with  all  GCBA  algorithms. 

CPU  time:  Vertical  particle  transport 

The  normalized  CPU  times  for  the  vertical  particle  transport  problems  are  displayed  in  Fig. 
11.  As  in  the  horizontal  case,  the  CPU  time  increases  with  increasing  grid  density.  For  the 
gas-solid  flow  problems  (Figs.  11(a)  and  11(b)),  the  relative  performance  of  the  various 
algorithms  is  similar  for  both  dilute  and  dense  concentration  of  particles.  For  the  dilute  case 
(Fig.  11(a)),  the  efficiency  of  PRIME  is  slightly  better  than  SIMPLEST  (due  to  the  use  of  an 
explicit  algebraic-equation  solver),  both  however  are  about  four  times  more  expensive  than 
all  other  algorithms  whose  performance  is  very  comparable  (i.e.  of  the  same  order  of 
magnitude)  with  SIMPLEM  being  the  least  expensive  (7%  less  than  SIMPLE  on  the  finest 
mesh)  and  PISO  the  most  expensive  (9.5%  more  than  SIMPLE  on  the  finest  mesh).  The  same 
is  true  for  the  dense  gas-solid  case  (Fig.  11(b))  with  the  performance  of  SIMPLE,  SIMPLEC, 
SIMPLEM,  SIMPLEX,  and  PISO  being  closer. 

For  the  vertical  bubbly  flows,  a  noticeable  change  in  the  normalized  time  chart  (Figs.  11(c) 
and  11(d))  is  depicted,  with  the  performance  of  SIMPLEST  and  PRIME  showing  good 
improvements  while  the  performance  of  the  remaining  algorithms  deterioration.  As  depicted, 
the  CPU  times  needed  by  the  various  algorithms  are  of  the  same  order  of  magnitude  with 
SIMPLEM  being  slightly  more  expensive. 
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By  comparing  the  behavior  of  the  various  algorithms  in  all  problems,  it  is  clear  that  the 
performance  of  SIMPLE,  SIMPLEC,  SIMPLEM,  SIMPLEX,  and  PISO  is  consistent  and 
require,  on  average,  the  least  computational  effort.  The  performance  of  the  SIMPLEST  and 
PRIME  algorithms  was  comparable  to  SIMPLE  for  upward  bubbly  flows  only  and  they  were, 
in  general,  the  most  expensive  to  use  on  all  grids  and  for  all  physical  situations  presented 
here.  Most  importantly  however,  is  the  fact  that  all  these  algorithms  can  be  used  to  predict 
multi-phase  (in  this  case  two-phase)  flows. 

Two-dimensional  two-phase  validation  problems 

In  this  section,  four  two-dimensional  two-phase  flow  problems  are  solved.  The  first  two 
problems  deal  with  incompressible  turbulent  flows  while  the  last  two  problems  are  concerned 
with  compressible  flows. 

Problem  1  .-Turbulent  upward  bubbly  flow  in  a  pipe 

Many  experimental  and  numerical  studies  involving  the  prediction  of  radial  phase 
distribution  in  turbulent  upward  air- water  flow  in  a  pipe  have  appeared  in  the  literature  [61- 
68].  These  studies  indicated  that  the  lateral  forces  that  most  strongly  affect  the  void 
distribution  are  the  lateral  lift  force  and  the  turbulent  stresses.  As  such,  in  addition  to  the 
usual  drag  force,  the  lift  force  is  considered  as  part  of  the  interfacial  force  terms  in  the 
momentum  equations.  In  the  present  work,  the  interfacial  drag  forces  per  unit  volume  are 
given  by: 

(ImE’  =-(lX=0375%'VdVW,lip(u'd’-u<->)  (132) 

r 

p 

(PS  =-(lyX=0375^p‘VdVdXlp(v<d'-v<d>)  (133) 

P 

where  rp  is  the  bubble  radius.  The  drag  coefficient  Cb  varies  as  a  function  of  the  bubble 
Reynolds  and  Weber  numbers  defined  as: 
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where  G,  the  surface  tension,  is  assigned  a  value  of  0.072  N/m  for  air-water  systems.  The 
drag  coefficient  is  computed  using  the  following  correlations  [69,70]: 


cD  =  - 


for  Re  „  <  0.49 


D  Re  p 643 

6.3 

D  Re0385 


for  0.49  <  Re„  <100 


for  Re  „  »  100 


(135) 


for  Re  p  » 100  and  We  >  8 


for  Re  p  »  100  and  Rep  >  2065.1/  We * 


Many  investigators  have  considered  the  modeling  of  lift  forces  [69-71].  Based  on  their  work, 
the  following  expressions  are  employed  for  the  calculation  of  the  interfacial  lift  forces  per 
unit  volume: 


GJLC>  =-(l„r  =ClP«»r«>(u»  -u">)x(VxaM) 
where  Ci  is  the  interfacial  lift  coefficient  calculated  from: 


(136) 


C,  =Cla(l-  2.78(0.2,  r(d))) 

where  (a,  b) denotes  the  minimum  of  a  and  b  and  Cia  is  an  empirical  constant. 


(137) 


The  effect  of  bubbles  on  the  turbulent  field  is  very  important.  In  this  work,  turbulence  is 
assumed  to  be  a  property  of  the  continuous  liquid  phase  (c)  and  is  computed  by  solving  Eqs. 
(4)  and  (5)  with  I®  and  1^’ given  by  [69]: 


I-  =V.  p(c)  3^  k(c,Vr(c)  +r(c)Pb 

o. 


(  A,(C)  \ 

!lc)=V.  p(c)  e(c,Vr(c)  +r(c)cl£Ph 


(138) 
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where  Gr  is  the  turbulent  Schmidt  number  for  volume  fractions,  and  Pb  is  the  production  rate 
of  k(c)  by  drag  due  to  the  motion  of  the  bubbles  through  the  liquid  and  is  given  by: 


Pb  = 


0.375CbCDp(c,r(d)r(c,V. 


2 

slip 


(140) 


In  Eq.  (140)  Q,  is  an  empirical  constant  representing  the  fraction  of  turbulence  induced  by 
bubbles  that  goes  into  large-scale  turbulence  of  the  liquid  phase.  Moreover,  as  suggested  in 
[69],  the  flux  representing  the  interaction  between  the  fluctuating  velocity  and  volume 
fraction  is  modeled  via  a  gradient  diffusion  approximation  and  added  as  a  source  term  in  the 

continuity  (v.(plk,D(k)Vr(k)))  and  momentum  (v.(p(k)D(k,u<kVr(k)))  equations  with  the 
diffusion  coefficient  D  given  by: 

v <k) 

D(k)  =  — —  (141) 


The  turbulent  viscosity  of  the  dispersed  air  phase  (d)  is  related  to  that  of  the  continuous  phase 
through: 


v  <c-* 

Vt(d)  =-L-  (142) 

Of 

where  of  is  the  turbulent  Schmidt  number  for  the  interaction  between  the  two  fluids.  The 
above  described  turbulence  model  is  a  modified  version  of  the  one  described  in  [69]  in  which 
the  turbulent  viscosities  of  both  fluids  are  allowed  to  be  different  in  contrast  to  what  is  done 
in  [69].  This  is  accomplished  through  the  introduction  of  the  of  parameter.  As  such,  different 
diffusion  coefficients  (D(k))  are  used  for  the  different  fluids.  Results  are  compared  against  the 
experimental  data  reported  by  Seriwaza  et  al  [61]. 

In  the  Seriwaza  et  al  experiment  [61],  the  Reynolds  number  based  on  superficial  liquid 
velocity  and  pipe  diameter  is  8xl04,  the  inlet  superficial  gas  and  liquid  velocities  are  0.077 
and  1.36  m/s,  respectively,  and  the  inlet  void  fraction  is  5.36x10  “  with  no  slip  between  the 
incoming  fluids.  Moreover,  the  bubble  diameter  is  taken  as  3  mm  [69],  while  the  fluid 
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properties  are  taken  as  p(c)=1000  Kg/m3,  p(d)=1.23  Kg/m3,  and  v[c)=  10'6  m2/s.  The  constants 
in  the  model  were  set  to:  Cia=0.075,  of=0.5,  or=0.7,  and  Cb=0.05.  Predicted  radial  profiles  of 
the  vertical  liquid  velocity  and  void  fraction  presented  in  Fig.  12  using  a  grid  of  size  96x32 
control  volumes  concur  very  well  with  measurements  and  compare  favorably  with  numerical 
profiles  reported  by  Boisson  and  Malin  [69].  As  shown,  the  void  fraction  profile  indicates 
that  gas  is  taken  away  from  the  pipe  center  towards  the  wall.  This  is  caused  by  the  lift  force, 
which  drives  the  bubbles  towards  the  wall. 

Having  established  the  credibility  of  the  physical  model  and  numerical  procedure,  the  next 
task  is  to  assess  the  merits  of  the  various  algorithms  for  such  flows.  For  that  purpose 
calculations  are  performed  using  the  SG,  PG,  and  FMG  strategies  for  all  algorithms.  Results 
are  displayed  in  the  form  of  (i)  total  mass  residuals  summed  over  both  phases  as  a  function  of 
outer  iterations  (Fig.  (13)),  and  (ii)  normalized  CPU  time  (Table  1)  needed  for  the  maximum 
normalized  residual  of  all  variables  and  for  all  phases  to  drop  below  £s=10~6. 

As  can  be  seen  from  Fig.  13,  it  is  possible  to  converge  the  solution  to  the  desired  level  with 
all  algorithms.  With  the  exception  of  PISO  (Fig.  13(a)),  the  convergence  characteristics  of  all 
algorithms  (Figs.  13(b)-(g))  are  very  similar  with  PRIME,  as  expected,  requiring  the  largest 
number  of  iterations.  The  PG  method  reduces,  on  average,  the  number  of  iterations  in 
comparison  with  the  SG  method  by  about  20%.  On  the  other  hand,  the  FMG  method  results 
in  a  50%  reduction  in  the  number  of  outer  iterations.  The  use  of  3  and  4  levels  for  both  PG 
and  FMG  methods  does  not  seem  to  have  any  effect  on  solution  acceleration  for  all 
algorithms  except  PISO,  for  which  the  use  of  4  levels  with  the  FMG  method  increases  the 
number  of  outer  iterations  considerably  and  results  in  a  kind  of  oscillations  (Fig.  13(a)).  The 
convergence  histories  of  all  algorithms  with  the  FMG  method  on  3  levels  presented  in  Fig. 
13(h)  confirm  once  more  the  aforementioned  observations. 


me  geometric  u unservauon  nasea  /ugoritnms  jor  iviuui-r  luia  now  at  /ut  jpeeas 


“-r-t 


In  Table  1,  the  normalized  CPU-times  (i.e.  CPU-time  divided  by  the  time  needed  by  SIMPLE 
on  the  coarsest  grid)  required  by  the  different  algorithms  to  converge  the  solution  to  the 
desired  level  with  the  various  methodologies  are  presented.  For  the  SG  method,  the  CPU¬ 
time  on  two  different  grids  of  sizes  48x16  and  96x32  C.V.  are  presented.  As  expected,  the 
CPU  effort  increases  for  all  algorithms  with  increasing  the  grid  size.  The  PG  and  FMG 
solutions  are  for  a  grid  of  96x32  C.V.  using  3  (96x32,  48x16,  and  24x8  C.V.)  and  4  (96x32, 
48x16,  24x8,  and  12x4  C.V.)  grid  levels. 

The  CPU  times  of  SIMPFE,  SIMPFEC,  and  SIMPFEX  are  very  close  on  all  grids  and  for  all 
methods  with  no  clear  superiority  of  any  algorithm  over  the  others.  The  performance  of 
SIMPFEST,  SIMPEEM,  and  PRIME  is  very  close  with  SIMPLEST  and  PRIME  being  the 
least  expensive  using  the  SG  and  PG  methods,  respectively.  With  the  FMG  method,  the 
PRIME  algorithm  is  the  most  expensive  with  no  clear  superiority  of  SIMPLEST  over 
SIMPLEM  and  vice  versa.  The  performance  of  PISO  with  the  SG  and  PG  methods  is  very 
close  to  that  of  the  SIMPLE  algorithm.  With  the  4  levels  FMG  method  however,  its 
performance  is  highly  unexpected  necessitating  higher  computational  effort  than  the  SG 
method  and  may  be  caused  by  the  additional  explicitness  introduced  by  the  PRIME  step.  The 
use  of  the  multi-grid  method  reduces  the  computational  cost,  on  average,  by  about  45%  (i.e. 
the  FMG  method  is  almost  twice  as  fast  as  the  SG  method)  while  the  use  of  the  PG  method 
results  in  a  reduction  of  about  15%,  both  in  comparison  with  the  SG  method.  Moreover,  with 
the  FMG  method,  the  least  computational  effort  is  obtained  with  SIMPLEC,  which  is  about 
0.7%  less  expensive  than  SIMPLE.  Excluding  PISO,  the  most  expensive  algorithm  with  the 
FMG  is  PRIME,  which  requires  27%  more  time  than  SIMPLE. 

Problem  2  .-Turbulent  air-particle  flow  in  a  vertical  pipe 

Here,  the  upward  flow  of  a  dilute  gas-solid  mixture  in  a  vertical  pipe  is  simulated.  As  in  the 
previous  problem,  the  axi- symmetric  form  of  the  gas  and  particulate  transport  equations  are 
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employed.  As  reported  in  several  studies  [72-74],  the  effects  of  interfacial  virtual  mass  and 
lift  forces  are  small  and  may  be  neglected  and  the  controlling  interfacial  force  is  drag 
(Harlow  and  Amsden  [75]),  which  is  given  by: 
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where  rp  represents  the  particle’s  radius,  Cd  the  drag  coefficient  computed  from: 
24 
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CD  =  0.44 
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and  Rep  the  Reynolds  number  based  on  the  particle  size  as  defined  in  Eq.(134). 

As  before,  turbulence  is  assumed  to  be  a  property  of  the  continuous  gas  phase  (c)  and  is 
predicted  using  a  two-fluid  k-e  model.  Several  extensions  of  the  k-8  model  for  carrier-phase 
turbulence  modulation  have  been  proposed  in  the  literature  [50-55]  and  the  one  suggested  by 
Chen  and  Wood  [52],  which  introduces  additional  source  terms  into  the  turbulence  transport 
equations,  is  adopted  here.  Thus,  the  turbulent  viscosity  is  computed  by  solving  the 
turbulence  transport  equations  (Eqs.  (4)  and  (5))  for  the  continuous  phase  with 
lj,k)  and  ((‘evaluated  using  the  following  relations  [52]: 
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where  xp  and  xe  are  timescales  characterizing  the  particle  response  and  large-scale  turbulent 
motion,  respectively,  and  are  computed  from: 
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with  Fd  being  the  magnitude  of  the  inter-fluid  drag  force  per  unit  volume.  The  turbulent  eddy 
viscosity  of  the  dispersed  phase  (d)  is  considered  to  be  a  function  of  that  of  the  continuous 
phase  and  is  computed  using  Eq.  (142). 

The  above-described  model  is  validated  against  the  experimental  results  of  Tsuji  et  al  [72]. 
Results  are  replicated  here  for  the  case  of  an  air  Reynolds  number,  based  on  the  pipe  diameter 
(of  value  30.5  mm),  of  3.3xl04  and  a  mean  air  inlet  velocity  of  15.6  m/s  using  particles  of 
diameter  200  pm  and  density  1020  Kg/m3.  In  the  computations,  the  mass-loading  ratio  at 
inlet  is  considered  to  be  1  with  no  slip  between  the  fluids,  and  Of  and  or  are  set  to  5  and  1010, 
respectively  (i.e.  the  interaction  terms  included  for  bubbly  flows  are  neglected  here).  Figure 
14  shows  the  fully  developed  gas  and  particles  mean  axial  velocity  profiles  generated  using  a 
grid  of  size  96x40  C.V.  It  is  evident  that  there  is  generally  a  very  good  agreement  between 
the  predicted  and  experimental  data  with  the  gas  velocity  being  slightly  over  predicted  and 
the  particles  velocity  slightly  under  predicted.  Moreover,  close  to  the  wall,  the  model 
predictions  indicate  that  the  particles  have  higher  velocities  than  the  gas,  which  is  in  accord 
with  the  experimental  results  of  Tsuji  et  al.  [72]. 

Having  checked  the  correctness  of  the  physical  model  and  numerical  procedure,  the  problem 
is  solved  using  the  SIMPLE,  SIMPLEC,  SIMPLEX,  and  SIMPLEST  multi-fluid  algorithms 
and  the  SG,  PG,  and  FMG  solution  methods.  As  in  the  previous  problem,  results  are 
displayed  in  the  form  of  (i)  total  mass  residuals  summed  over  both  phases  as  a  function  of 
outer  iterations  (Fig.  (15)),  and  (ii)  normalized  CPU  time  (Table  2)  needed  for  the  maximum 
normalized  residual  of  all  variables  and  for  all  phases  to  drop  below  £s=10~6. 
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Mass  residual  plots  presented  in  Fig.  15  indicate  a  similar  convergence  behavior  for  all  four 
algorithms  with  very  close  number  of  outer  iterations  to  achieve  the  desired  level  of 
convergence  (i.e.  within  50  outer  iterations).  It  is  hard  to  see  any  noticeable  difference 
between  the  3  and  4  levels  with  both  the  PG  and  FMG  methods.  The  decrease  in  the  number 
of  iterations  with  the  PG  method  over  the  SG  method  is  smaller  than  the  decrease  obtained 
with  bubbly  flows.  This  lower  effectiveness  of  the  PG  method  is  due  to  the  following  reason. 
In  solving  the  problem,  it  is  noticed  that  the  initial  guess  greatly  affects  the  convergence 
history  and  time  required  to  converge  the  solution  to  the  desired  level.  Except  when  solving 
on  the  finest  mesh  with  the  SG  method,  the  initial  guess  used  for  the  velocity  field  is 
u(c)=u<d,=l  m/s.  The  use  of  this  initial  value  with  the  SG  method  on  the  finest  mesh  greatly 
increased  the  CPU  effort  needed  over  the  one  needed  when  starting  with  an  initial  field  of 
u(c)=u<d,=15.6  m/s.  To  reduce  cost,  the  latter  initial  guess  is  used.  For  this  reason  the  mass 
residuals  start  from  somehow  a  lower  value  than  expected  and  the  PG  method  appears  to  be 
less  effective.  The  FMG  method  reduces  the  number  of  outer  iterations  by  about  53%  over 
the  SG  method,  which  indicates  a  good  capability  to  deal  with  the  added  non-linearity  of 
multiphase  flows. 

The  normalized  CPU-times  presented  in  Table  2  reflects  the  above  stated  behaviors  with  the 
time  required  when  using  the  PG  method  being  very  close  (slightly  higher,  except  with 
SIMPLEM  it  is  slightly  lower)  to  the  time  needed  with  the  SG  method.  The  FMG  method 
reduces  the  cost,  on  average,  by  about  40  %.  The  least  computational  effort  is  accomplished 
with  SIMPLE  while  SIMPLEM  is  the  most  expensive  with  all  methods  (23%  more  expensive 
than  SIMPLE  with  the  FMG  on  3  levels). 

Problem  3 .‘Compressible  dilute  air-particle  flow  over  a  flat  plate 

As  has  been  demonstrated  in  several  studies  [76-82],  two-fluid  flow  greatly  changes  the  main 

features  of  the  boundary  layer  over  a  flat  plate.  Typically,  three  distinct  regions  are  defined  in 
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the  two-fluid  boundary  layer  (Fig.  16),  based  on  the  importance  of  the  slip  velocity  between 
the  two  fluids:  a  large-slip  region  close  to  the  leading  edge,  a  moderate-slip  region  further 
down,  and  a  small-slip  zone  far  downstream.  The  characteristic  scale  in  this  two-fluid  flow 
problem  is  the  relaxation  length  A,e  [79],  defined  as: 


V  = 


2  P<d)rp2uc 


(149) 


9  p(c) 

where  u  ^  is  the  free  stream  velocity.  The  three  regions  are  defined  according  to  the  order  of 
magnitude  of  the  slip  parameter  x*  =  x/  \ .  In  the  simulation,  the  viscosity  of  the  fluid  is 
considered  to  be  a  function  of  temperature  according  to  [79]: 


(  rp<c>  \°'6 


(150) 


V  "ref  J 


where  the  reference  viscosity  and  temperature  are  pref=1.86xl0  5  N.s/m2  and  Tref=303  °K. 
Drag  is  the  only  retained  interfacial  force  as  it  dominates  the  other  interfacial  forces.  It  is 
computed  as  [79]: 
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where  the  drag  coefficient  is  given  by: 

CD  =—  Rep  +  -Re°15 
50  P  6  P 


(151) 
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(153) 


In  the  energy  equation,  heat  transfer  due  to  radiation  is  neglected  and  only  convective  heat 
transfer  around  an  isolated  particle  is  considered.  Under  such  conditions,  the  interfacial  terms 
in  the  gas  (c)  and  particles  (d)  energy  equations  reduce  to  [79]: 


I 


(C) 

E 

l.d> 


Qg-p +Fg-p-u 
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r  =  -o 

LE  Vg-p 


(154) 
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where: 
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In  the  above  equations,  Nu  is  the  Nusselt  number,  Pr(c)  the  gas  Prandtl  number,  X<c)  the  gas 
thermal  conductivity,  T  the  temperature,  and  other  parameters  are  as  defined  earlier. 

In  the  simulation,  the  particle  diameter,  particle  Reynolds  number,  material  density,  Prandtl 
number,  and  mass  load  ratio  are  set  to:  10  (im,  10,  1766  kg/m3,  0.75,  and  1  respectively.  The 
wall  boundary  is  treated  as  a  no-slip  boundary  for  the  gas  phase  (i.e.  both  components  of  the 
gas  velocity  are  set  to  zero),  and  as  a  slip  boundary  condition  for  the  particles  phase  (i.e.  the 
normal  fluxes  are  set  to  zero).  In  order  to  bring  all  quantities  to  the  same  order  of  magnitude, 
results  are  displayed  using  the  following  dimensionless  variables: 


x*  =  —  ,y*  =  -^VR^,u*=  — ,v*  =  —  VRe  Re  =  ^  (159) 

Figure  17  shows  the  results  for  the  steady  flow  obtained  on  a  rectangular  domain  with  a  mesh 
of  density  104x48  C.V.  stretched  in  the  y-direction.  The  figure  depicts  the  development  of 
gas  and  particles  velocity  profiles  within  the  three  regions  mentioned  earlier.  In  the  near 
leading  edge  area  (x*=0.1),  the  gas  velocity  is  adjusted  at  the  wall  to  obtain  the  no-slip 
condition  as  for  the  case  of  a  pure  gas  boundary  layer.  The  particles  have  no  time  to  adjust  to 
the  local  gas  motion  and  there  is  a  large  velocity  slip  between  the  fluids.  In  the  transition 
region  (x*=l),  significant  changes  in  the  flow  properties  take  place.  The  interaction  between 
the  fluids  cause  the  particles  to  slow  down  while  the  gas  accelerates.  In  the  far  downstream 
region  (x*=5),  the  particles  have  ample  time  to  adjust  to  the  state  of  the  gas  motion,  the  slip  is 
very  small,  and  the  solution  tends  to  equilibrium.  These  results  are  in  excellent  agreement 
with  numerical  solutions  reported  by  Thevand  et  al.  [82]  (Fig.  17). 
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As  in  test  1,  the  problem  is  solved  using  all  multi-fluid  algorithms  and  the  SG,  PG,  and  FMG 
solution  methods.  Results  are  displayed  in  the  form  of  total  mass  residuals  (Fig.  18)  and 
normalized  CPU  time  (Table  3)  with  £s=10"6. 

Plots  presented  in  Fig.  18  indicate  that  it  is  possible  to  converge  the  solution  to  the  desired 
level  with  all  algorithms.  As  depicted  in  figure  18(a),  PISO  requires  the  least  number  of  outer 
iterations.  This,  however,  is  not  associated  with  the  lowest  computational  effort  due  to  the 
higher  cost  per  iteration  in  comparison  with  other  algorithms.  The  convergence 
characteristics  of  SIMPLE  (Fig.  18(b)),  SIMPLEC  (Fig.  18(c)),  SIMPLEM  (Fig.  18(d)),  and 
SIMPLEX  (Fig.  18(g))  are  very  similar,  requiring  nearly  the  same  number  of  outer  iterations 
with  the  SG,  PG,  and  FMG  methods.  The  number  of  iterations  required  by  PRIME  (Fig. 
18(f))  with  the  SG  method  is  higher  than  SIMPLEST  (Fig.  18(e)).  With  the  FMG  method 
however,  the  performance  of  the  two  algorithms  is  very  close  with  that  of  PRIME  being 
slightly  better.  In  general,  the  use  of  the  PG  method  reduces  the  number  of  outer  iterations,  as 
compared  to  the  SG  method,  by  over  40%  with  all  algorithms  whereas  the  use  of  the  FMG 
method  reduces  it  by  over  64%.  Fig.  18(h)  indicates  that  when  using  the  FMG  method  on  3 
levels  PISO  requires  the  lowest  number  of  iterations  followed  by  SIMPLEM,  SIMPLEC, 
SIMPLEX,  SIMPLE,  PRIME,  and  SIMPLEST.  With  the  number  of  iterations  needed  by 
SIMPLEST  and  PRIME  being  very  close  and  nearly  double  the  number  needed  by  SIMPLE. 
It  should  be  clarified  that  the  displayed  numbers  of  iterations  represent  those  needed  for  the 
mass  residuals  to  be  reduced  to  the  desired  level.  The  CPU-times  however  represent  the 
computational  effort  needed  to  reduce  the  maximum  residuals  to  the  desired  level.  In  some 
cases,  even  though  the  mass  residuals  may  become  below  the  desired  value,  other  residuals 
could  still  be  above  that  value.  This  is  why,  for  example,  the  CPU  time  needed  by  SIMPLE  is 
lower  than  that  needed  by  SIMPLEC  (Table  3)  even  though  Fig.  18(h)  indicates  that  the 
number  of  iterations  needed  by  SIMPLEC  to  reduce  the  mass  residuals  to  below  8S  is  lower. 
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In  Table  3,  the  normalized  CPU-times  required  by  the  different  algorithms  to  solve  the 
problem  with  the  various  methodologies  are  presented.  For  the  SG  method,  the  CPU-time  on 
two  different  grids  of  sizes  52x24  and  104x48  C.V.  are  presented.  As  expected,  the  CPU 
effort  increases  for  all  algorithms  with  increasing  the  grid  size.  The  PG  and  FMG  solutions 
are  for  a  grid  of  104x48  C.V.  using  3  (104x48,  52x24,  and  26x12  C.V.)  and  4  (104x48, 
52x24,  26x12,  and  13x6  C.V.)  grid  levels.  The  SIMPLE  algorithm  appears  to  be  the  most 
efficient  requiring  the  least  computational  effort  on  all  grids  and  with  all  methods.  The 
performance  of  SIMPLEC  and  SIMPLEX  is  the  closest  to  SIMPLE  especially  with  the  multi¬ 
grid  method.  SIMPLEST  and  PRIME  are  the  most  expensive  to  use  with  SIMPLEST 
performing  better  with  the  SG  and  PG  methods  (PRIME  requires  170%  more  time  than 
SIMPLE  on  the  dense  grid  with  the  SG  method  while  SIMPLEST  requires  54%  additional 
time).  PISO  is  more  expensive  than  SIMPLE  on  all  grids  and  with  all  methods  (it  requires 
15%  more  time  than  SIMPLE  on  the  dense  grid  with  the  SG  method  and  56%  more  time  with 
the  FMG  on  3  levels). 

Problem  4:  Inviscid  transonic  dusty  flow  in  a  converging-diverging  nozzle 

The  last  test  considered  deals  with  the  prediction  of  supersonic  dilute  air-particle  flow  in  an 

axi-symmetric  converging-diverging  rocket  nozzle.  Several  researchers  have  analyzed  the 
problem  and  data  is  available  for  comparison  [83-92].  In  most  of  the  reported  studies,  a 
shorter  diverging  section,  in  comparison  with  the  one  considered  here,  has  been  used  when 
predicting  the  two-fluid  flow.  Two-fluid  flow  results  for  the  long  configuration  have  only 
been  reported  by  Chang  et.al.  [87].  The  flow  is  assumed  to  be  inviscid  and  the  single-fluid 
flow  results  are  used  as  an  initial  guess  for  solving  the  two-fluid  flow  problem.  The  physical 
configuration  (Fig.  19)  is  the  one  described  in  [87].  The  viscosity  of  the  fluid  varies  with  the 
temperature  according  to  Sutherland’s  law  for  air: 
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The  coupling  between  gas  and  particle  phases  is  through  the  interfacial  momentum  and 
energy  terms.  The  force  exerted  on  a  single  particle  moving  through  a  gas  [88]  is  given  as: 

fx  =6ra-pfD(ilc)(u<d) -u<c))  (161) 

fy  =67trpfD|i(c)(v(d)-v(c))  (162) 

so  that  for  N  particles  in  a  unit  volume  the  effective  drag  force  is 


(r  )(c)  =  _(r  )' 

V^m  / d  \  m  )\ 


l(d)  9  r(d) 


2  r 


fDH*>(uM>-uM) 


(rf =-(i’X  -vw) 


p 

(d) 


2  r ' 


(163) 


(164) 


where  to  is  the  ratio  of  the  drag  coefficient  Cd  to  the  stokes  drag  CDo=24/Rep  and  is  given  by 
[87] 


fD  =  1  +  0.15Rep687+ 


0.0175Re 

- - — - — — —  Re  <  3x10 

1  +  4.25xl04  Re“  16  p 


(165) 


The  heat  transferred  from  gas  to  particle  phase  per  unit  volume  is  given  as  [88] 


3  r(d) 

Q  =- — 

s'p  2  rD 


X(c)Nu(T(d)  -  T(c) ) 


(166) 


Where  X(c)  is  the  thermal  conductivity  of  the  gas  and  Nu  the  Nusselt  number,  which  is  written 
as  [88] 


Nu  -2  +  0.459 Re “ 55  Pr"33  (167) 

The  gas-particle  inter-fluid  energy  term  is  given  by 

1“  =  ~fDn"’(u|d|  -uw)ua  +~fDM«(vla>  -  v«)va  +  X.ftlNu(T|d>-TM)  (168) 
2  rp  2  rp  2  rp 

I<d)  =  l£lAx(c>Nu(T(0)  -  T(d) )  (169) 

2  rp 

where  the  first  two  terms  on  the  right-hand  side  of  equation  (168)  represent  the  energy 


exchange  due  to  momentum  transfer. 
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The  physical  quantities  employed  are  similar  to  those  used  in  [87].  The  gas  stagnation 
temperature  and  pressure  at  inlet  to  the  nozzle  are  555  °K  and  10.34x10'  N/m  ,  respectively. 
The  specific  heat  for  the  gas  and  particles  are  1.07xl03  J/Kg°K  and  1.38xl03  J/Kg°K, 
respectively,  and  the  particle  density  is  4004.62  kg/m  .  With  a  zero  inflow  velocity  angle,  the 
fluid  is  accelerated  from  subsonic  to  supersonic  speed  in  the  nozzle.  The  inlet  velocity  and 
temperature  of  the  particles  are  presumed  to  be  the  same  as  those  of  the  gas  phase.  Results  for 
two  particle  sizes  of  radii  1  and  10  (im  with  the  same  mass  fraction  cf)=0.3  are  presented  using 
a  grid  of  size  188x80  C.V.  Figures  20(a)  and  20(b)  show  the  particle  volume  fraction 
contours  while  Figures  20(c)  and  20(d)  display  the  velocity  distribution.  For  the  flow  with 
particles  of  radius  1  (im,  a  sharp  change  in  particle  density  is  obtained  near  the  upper  wall 
downstream  of  the  throat,  and  the  particle  density  decreases  to  a  small  value.  With  the  large 
particle  flow  (10  (im),  however,  a  much  larger  particle-free  zone  appears  due  to  the  inability 
of  the  heavier  particles  to  turn  around  the  throat  corner.  These  findings  are  in  excellent 
agreement  with  published  results  reported  in  [87]  and  others  using  different  methodologies.  A 
quantitative  comparison  of  current  predictions  with  published  experimental  and  numerical 
data  is  presented  in  Fig.  21  through  gas  Mach  number  distributions  along  the  wall  (Fig.  21(a)) 
and  centerline  (Fig.  21(b))  of  the  nozzle  for  the  one-fluid  and  two-fluid  flow  situations.  As 
can  be  seen,  the  one-fluid  flow  predictions  fall  on  top  of  experimental  data  reported  in  [90- 
92].  Since  the  nozzle  contour  has  a  rapid  contraction  followed  by  a  throat  with  a  small  radius 
of  curvature,  the  flow  near  the  throat  wall  is  overturned  and  inclined  to  the  downstream  wall. 
A  weak  shock  is  thus  formed  to  turn  the  flow  parallel  to  the  wall.  This  results  in  a  sudden 
drop  in  the  Mach  number  value  and  as  depicted  in  Fig.  21(b),  this  sudden  drop  is  correctly 
envisaged  by  the  solution  algorithm  with  the  value  after  the  shock  being  slightly  over 


predicted. 
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Due  to  the  unavailability  of  experimental  data,  two-fluid  flow  predictions  are  compared 
against  the  numerical  results  reported  in  [87].  As  displayed  in  Figs.  21(a)  and  21(b),  both 
solutions  are  in  good  agreement  with  each  other  indicating  once  more  the  correctness  of  the 
calculation  procedures.  The  lower  gas  Mach  number  values  in  the  two-fluid  flow  is  caused  by 
the  heavier  particles  <jp(d)»p(c)),  which  reduce  the  gas  velocity.  Moreover,  owing  to  the 
particle-free  zone,  the  Mach  number  difference  between  the  one-  and  two-fluid  flows  along 
the  wall  is  smaller  than  that  at  the  centerline. 

To  compare  the  relative  performance  of  the  multi-fluid  algorithms,  the  problem  is  solved  via 
the  PG  method  using  the  SIMPLE,  SIMPLEC,  SIMPLEM,  and  SIMPLEX  algorithms  over 
three  different  grids  of  sizes  47x20,  94x40,  and  188x80  C.V.  for  a  particle  radius  of  size  1 
pm.  As  before,  results  are  displayed  in  the  form  of  total  mass  residuals  (Fig.  22)  and 
normalized  CPU  times  (Table  4)  with  £s  set  to  10"5.  As  shown  in  Fig.  22,  all  algorithms 
require  almost  the  same  number  of  iterations  with  the  exception  of  SIMPLE  on  the  94x40 
grid,  which  requires  a  larger  number  of  iterations  than  on  the  finest  mesh.  Excluding  that 
case,  the  convergence  histories  of  all  algorithms  are  nearly  identical.  In  terms  of 
computational  effort,  the  normalized  CPU-times  presented  in  Table  4  indicate  a  close 
performance  of  the  various  algorithms  with  SIMPLE  being  the  most  efficient  (7%  less 
expensive  than  SIMPLEC)  and  SIMPLEM  the  most  expensive  (43%  more  expensive  than 
SIMPLE)  on  the  finest  mesh.  On  the  other  hand,  SIMPLEX  is  11%  more  expensive  than 
SIMPLE. 

Closing  Remarks 

The  implementation  of  seven  GCBA  multi-fluid  algorithms  for  the  simulation  of  multi-fluid 
flow  at  all  speeds  was  accomplished.  The  algorithms  were  embedded  within  a  non-linear  full 
multi-grid  strategy.  A  two-fluid  k-£  model  and  several  inter-phase  models  were  also 
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employed.  Solving  a  variety  of  one-  and  two-dimensional  two-phase  flow  problems  assessed 
the  performance  and  accuracy  of  these  algorithms.  For  each  test  problem,  solutions  were 
generated  on  a  number  of  grid  systems  using  the  single  grid  method  (SG),  the  prolongation 
grid  method  (PG),  and  the  full  non-linear  multi-grid  method  (FMG).  Results  obtained 
demonstrated  the  capability  of  all  algorithms  to  deal  with  multi-fluid  flow  situations  and  to 
predict  multi-fluid  flow  at  all  speeds,  and  the  ability  of  the  FMG  method  to  tackle  the  added 
non-linearity  of  laminar  and  turbulent  multi-fluid  flows.  The  convergence  history  plots  and 
CPU-times  presented,  indicated  similar  performances  for  SIMPLE,  SIMPLEC,  and 
SIMPLEX.  The  PISO,  SIMPLEM,  and  SIMPLEST  algorithms  were  in  general  more 
expensive  than  SIMPLE.  In  general,  the  PRIME  algorithm  was  the  most  expensive  to  use. 
The  PG  and  FMG  methods  accelerated  the  convergence  rate  for  all  algorithms.  The  FMG 
method  was  found  to  be  more  efficient. 
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Figure  Captions 


Fig.  1  (a)  Control  volume,  (b)  the  prolongation  only,  and  (c)  FMG  strategies,  and 
(d)  Physical  domain  for  the  gas-particle  transport  problem. 

Fig.  2  (a)  Comparison  between  the  analytical  and  numerical  particle  velocity  distributions, 

(b) -(g)  convergence  histories  on  the  different  grid  systems,  (h)  and  convergence 
histories  on  the  80  C.V.  grid  for  the  horizontal  dilute  gas-solid  flow  problem. 

Fig.  3  (a)  gas  and  particle  velocity  distributions,  (b)-(g)  convergence  histories  on  the 

different  grid  systems,  (h)  and  convergence  histories  of  the  various  algorithms  on  the 
80  C.V.  grid  for  the  horizontal  dense  gas-solid  flow  problem. 

Fig.  4  (a)  Liquid  and  gas  velocity  distributions,  (b)-(g)  convergence  histories  on  the  different 
grid  systems,  (h)  and  convergence  histories  of  the  various  algorithms  on  the  80  C.V. 
grid  for  the  horizontal  dilute  bubbly  flow  problem. 

Fig.  5  (a)  Liquid  and  gas  velocity  distributions,  (b)-(g)  convergence  histories  on  the  different 
grid  systems,  (h)  and  convergence  histories  on  the  80  C.V.  grid  for  the  horizontal 
dense  bubbly  flow  problem. 

Fig.  6  Normalized  CPU-times  for  the  horizontal  (a)  dilute  gas-solid,  (b)  dense  gas-solid, 

(c)  dilute  bubbly,  and  (d)  dense  bubbly  flow  problem. 

Fig.  7  (a)  gas  and  particle  velocity  distributions,  and  (b)-(h)  convergence  histories  on  the 
different  grid  systems  for  the  vertical  dilute  gas- solid  flow  problem. 

Fig.  8  (a)  gas  and  particle  velocity  distributions,  and  (b)-(h)  convergence  histories  on  the 
different  grid  systems  for  the  vertical  dense  gas-solid  flow  problem. 

Fig.  9  (a)  gas  and  particle  velocity  distributions;  and  (b)-(h)  convergence  histories  on  the 
different  grid  systems  for  the  vertical  dilute  bubbly  flow  problem. 

Fig.  10  (a)  gas  and  particle  velocity  distributions,  and  (b)-(h)  convergence  histories  on  the 
different  grid  systems  for  the  vertical  dense  bubbly  flow  problem. 

Fig.  1 1  Normalized  CPU-times  for  the  vertical  (a)  dilute  gas-solid,  (b)  dense  gas-solid, 


(c)  dilute  bubbly,  and  (d)  dense  bubbly  flow  problem. 

Fig.  12  Comparison  of  fully  developed  liquid  velocity  and  void  fraction  profiles  for  turbulent 
bubbly  upward  bubbly  flow  in  a  pipe. 

Fig.  13  (a)-(g)  Convergence  histories  of  the  SG,  PG,  and  FMG  methods  on  the  finest  grid, 
and  (h)  convergence  histories  of  the  various  algorithms  on  the  finest  mesh  using  the 
FMG  method  for  turbulent  upward  bubbly  flow  in  a  pipe. 

Fig.  14  Comparison  of  fully  developed  gas  and  particle  velocity  profiles  for  turbulent  air- 
article  flow  in  a  pipe. 

Fig.  15  Convergence  histories  of  the  (a)  SIMPLE,  (b)  SIMPLEC,  (c)  SIMPLEST,  and  (d) 
SIMPLEX  algorithms  using  the  SG,  PG,  and  FMG  methods  on  the  finest  mesh  for 
turbulent  air-particle  flow  in  a  pipe. 

Fig.  16  The  three  different  regions  within  the  boundary  layer  of  dusty  flow  over  a  flat  plate. 

Fig.  17  Comparison  of  fully  developed  gas  and  particle  velocity  profiles  inside  the  boundary 
layer  at  different  axial  locations  for  dilute  two-phase  flow  over  a  flat  plate. 

Fig.  18  (a)-(g)  Convergence  histories  of  the  SG,  PG,  and  FMG  methods  on  the  finest  grid, 
and  (h)  convergence  histories  of  the  various  algorithms  on  the  finest  mesh  using  the 
FMG  method  for  dusty  gas  flow  over  a  flat  plate. 

Fig.  19  Physical  domain  for  the  dusty  gas  flow  in  a  converging-diverging  nozzle. 

Fig.  20  (a,b)  Volume  Fraction  contours  and  (c,d)  particle  velocity  vectors  for  dusty  gas  flow 
in  a  converging-diverging  nozzle. 

Fig.  21  Comparison  of  one-phase  and  two-phase  gas  Mach  number  distributions  along  the  (a) 
wall  and  (b)  centerline  of  the  dusty  flow  in  a  converging-diverging  nozzle  problem. 

Fig.  22  Convergence  histories  of  the  (a)  SIMPLE,  (b)  SIMPLEC,  (c)  SIMPLEM,  and  (d) 
SIMPLEX  algorithms  using  the  SG  method  for  dusty  gas  flow  in  a  converging- 
diverging  nozzle. 
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Fig.  5  (a)  Liquid  and  gas  velocity  distributions,  (b)-(g)  convergence  histories  on  the  different  grid  systems,  (h)  and  convergence  histories  on  the  80 


C.V.  grid  for  the  horizontal  dense  bubbly  flow  problem. 
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ig.  7  (a)  gas  and  particle  velocity  distributions,  and  (b)-(h)  convergence  histories  on  the  different  grid  systems  for  the  vertical  dilute  gas-solid  flow 
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ig.  13  (a)-(g)  Convergence  histories  of  the  SG,  PG,  and  FMG  methods  on  the  finest  grid,  and  (h)  convergence  histories  of  the  various  algorithms  on  the 
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Fig.  14  Comparison  of  fully  developed  gas  and  particle  velocity  profiles  for  turbulent 


air-particle  flow  in  a  pipe. 
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Fig.  15  Convergence  histories  of  the  (a)  SIMPLE,  (b)  SIMPLEC,  (c)  SIMPLEM,  and 

(d)  SIMPLEX  algorithms  using  the  SG,  PG,  and  FMG  methods  on  the  finest  mesh  for 
turbulent  air-particle  flow  in  a  pipe. 
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Fig.  16  The  three  different  regions  within  the  boundary  layer  of  dusty  flow  over  a  flat 


plate. 


Fig.  17  Comparison  of  fully  developed  gas  and  particle  velocity  profiles  inside  the  boundary 


layer  at  different  axial  locations  for  dilute  two-phase  flow  over  a  flat  plate. 
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Fig.  18  (a)-(g)  Convergence  histories  of  the  SG,  PG,  and  FMG  methods  on  the  finest  grid,  and  (h)  convergence  histories  of  the  various  algorithms  on  the 


finest  mesh  using  the  FMG  method  for  dusty  gas  flow  over  a  flat  plate 
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Fig.  19  Physical  domain  for  the  dusty  gas  flow  in  a  converging-diverging  nozzle. 
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Fig.  20  (a,b)  Volume  Fraction  contours  and  (c,d)  particle  velocity  vectors  for  dusty  gas 
flow  in  a  converging-diverging  nozzle. 
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Fig.  21  Comparison  of  one-phase  and  two-phase  gas  Mach  number  distributions  along  the  (a) 


wall  and  (b)  centerline  of  the  dusty  flow  in  a  converging-diverging  nozzle  problem. 
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Fig.  22  Convergence  histories  of  the  (a)  SIMPLE,  (b)  SIMPLEC,  (c)  SIMPLEM,  and 


(d)  SIMPLEX  algorithms  using  the  PG  method  for  dusty  gas  flow  in  a  converging- 


diverging  nozzle. 


Table  1  Normalized  CPU-times  for  turbulent  bubbly  flow  in  a  pipe. 


ALGORITHMS 


GRID 

METHOD 

SIMPFE 

SIMPFEC 

SIMPFEX 

SIMPFEST 

SIMPEEM 

PISO 

PRIME 

48x16  C.V. 

SG 

1.00 

1.00 

1.02 

1.04 

1.09 

1.08 

1.19 

SG 

40.09 

40.42 

40.98 

40.66 

43.07 

41.82 

43.59 

PG  (3  levels) 

34.24 

34.41 

35.02 

36.82 

37.12 

34.48 

36.05 

96x32  C.V. 

PG  (4  levels) 

34.34 

34.25 

34.67 

36.88 

36.93 

34.50 

36.07 

FMG  (3  levels) 

22.32 

22.17 

22.40 

25.55 

23.67 

27.31 

28.38 

FMG  (4  levels) 

22.53 

22.51 

22.78 

23.81 

25.55 

55.68 

27.89 

Table  2  Normalized  CPU-times  for  turbulent  air-particle  flow  in  a  pipe. 

ALGORITHMS 


GRID 

METHOD 

SIMPFE 

SIMPFEC 

SIMPFEX 

SIMPEEM 

48x20  C.V. 

SG 

1.00 

1.26 

1.28 

1.45 

SG 

18.81 

18.85 

19.07 

21.05 

PG  (3  levels) 

18.94 

19.37 

19.37 

20.22 

96x40  C.V. 

PG  (4  levels) 

19.03 

19.46 

19.44 

20.30 

FMG  (3  levels) 

11.11 

12.34 

12.44 

13.72 

FMG  (4  levels) 

12.07 

13.10 

13.25 

13.54 

Table  3  Normalized  CPU-times  for  Dusty  flow  over  a  flat  plate. 


ALGORITHMS 


GRID 

METHOD 

SIMPFE 

SIMPFEC 

SIMPFEX 

SIMPFEST 

SIMPEEM 

PISO 

PRIME 

52x24  C.V. 

SG 

1.00 

1.04 

1.11 

1.52 

1.31 

1.12 

2.13 

SG 

16.92 

18.09 

18.94 

26.00 

21.48 

19.46 

45.87 

PG  (3  levels) 

14.25 

14.94 

15.51 

18.85 

17.48 

15.84 

24.48 

104x48  C.V. 

PG  (4  levels) 

14.27 

14.80 

15.50 

18.90 

17.56 

15.79 

24.46 

FMG  (3  levels) 

5.35 

5.44 

5.88 

11.87 

7.84 

8.34 

12.19 

FMG  (4  levels) 

6.12 

6.14 

6.66 

12.98 

6.66 

9.08 

9.98 

Table  4  Normalized  CPU-times  for  Dusty  flow  in  a  converging-diverging  nozzle. 

ALGORITHMS 


GRID 

METHOD 

SIMPLE 

SIMPLEC 

SIMPLEX 

SIMPLEM 

47x20  C.V. 

SG 

1.00 

1.02 

1.06 

1.13 

94x40  C.V. 

PG  (2  levels) 

13.08 

9.24 

9.47 

11.37 

188x80  C.V. 

PG  (3  levels) 

82.43 

88.49 

91.35 

117.83 

